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ABSTRACT 

Fractal  Geometry  is  a  recent  synthesis  of  old  mathematical  constructs.  It  was 
first  popularized  by  complex  renderings  of  terrain  on  a  computer  graphics 
medium.  Fractal  geometry  has  since  spawned  research  in  many  diverse  scientific 
disciplines.  Its  rapid  acceptance  has  been  achieved  due  to  its  ability  to  model 
phenomena  that  defy  discrete  computation  due  to  roughness  and  discontinuities. 
With  its  quick  acceptance  has  come  problems.  Fractal  geometry  is  a 
misunderstood  idea  that  is  quickly  becoming  buried  under  grandiose  terminology 
that  serves  no  purpose.  Its  essence  is  induction  using  simple  geometric  constructs 
to  transform  initiating  objects.  The  fractal  objects  that  we  create  with  this 
process  often  resemble  natural  phenomenon.  The  purpose  of  this  work  is  to 
present  fractal  geometry  to  the  graphics  programmer  as  a  simple  workable 
technique.  We  hope  to  demystify  the  concepts  of  fractal  geometry  and  make  it 
available  to  all  who  are  interested. 
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I.    AN  INTRODUCTION  TO  FRACTAL  GEOMETRY 

A.    MATHEMATICS  AS  A  MODEL  FOR  OUR  UNIVERSE 

The  various  branches  of  mathematics  have  through  time  developed  as  a 
response  to  the  need  for  more  detailed  models  to  describe  new  developments, 
both  technological  and  philosophical.  This  was  true  when  Newton  developed 
calculus  and  also  true  during  the  late  1800's  through  the  1920's  when  a  schism 
developed  between  the  classical  mathematicians  and  some  brilliant  innovative 
thinkers. 

1.     The  Mathematical  Crises  of  the  Early  19th  Century 

One  of  man's  greatest  strengths  is  his  ability  to  question  his  surroundings 
and  beliefs  and  through  this  questioning  develop  new  insight  and  innovation. 
Most  mathematical  systems  are  developed  for  use  in  applications.  Man's  natural 
inquisitiveness  often  leads  him  to  develop  his  systems  beyond  the  application  and 
into  abstract  theory.  This  theory  drives  him  to  investigate  the  applications  and 
often  yields  direction  for  new  discoveries  that  were  not  previously  foreseen  or  that 
defy  intuition. 

Georg  Cantor  (1845-1918)  was  the  most  notable  of  a  number  of 
mathematicians  who  questioned  the  basic  precepts  of  mathematics  and  developed 
the  modern  set  theory.  Some  of  Cantor's  discoveries  seemed  to  invalidate  many 
of  the  long  held  beliefs  of  mathematics.  Cantor  and  his  peers  became  deeply 
involved  in  controversy  over  their  findings.  Their  discovery  of  functions  which 
seemed  to  violate  the  basic  rules  of  geometry  and  calculus  were  deemed  as 
monsters  and  unworthy  of  consideration  by  reasonable  men  because  they  lacked 
usefulness  to  any  application  then  known  [Ref.  l:pp.  9].  These  new  concepts 
would  remain  in  the  arena  of  pure  theoretical  mathematics  until  science 
developed  to  a  point  where  the  old  models  could  no  longer  adequately  describe  it> 
processes  and  would  look  to  the  new  mathematics  for  a  new  perspective. 

It  was  from  these  discoveries  that  Fractal  Geometry  was  born  Ref. 
l:Chap.    2].     It    will    be    seen    in    the    following    chapters    that    fractal    geometry 


is  a  synthesis  of  many  of  the  concepts  which  developed  from  the  mathematical 
schism  of  the  19th  century,  most  notably  set  theory  and  topology. 

2.  What  is  a  Mathematical  Model 

Reference  2  defines  a  mathematical  model  in  the  following  fashion  [Ref. 
2:pp  1-3]: 

A  mathematical  model  is  a  mathematical  characterization  of  a  phenomenon  or 
process.  It  has  three  essential  parts:  a  process  or  phenomenon  which  is  to  be 
modeled,  a  mathematical  structure  capable  of  expressing  the  important  proper- 
ties of  the  object  to  be  modeled,  and  an  explicit  correspondence  between  the 
two. 

Although  the  phenomenon  of  interest  need  not  be  taken  from  the  real  world, 
they  usually  are.  The  real  world  component  is  described  quantitatively  by  such 
things  as  parameter  values  and  at  which  time  things  occur. 

The  second  component  of  a  model  is  an  abstract  mathematical  structure.  In  it- 
self, the  structure  is  abstract  and  has  no  intrinsic  relation  to  the  real  world. 
However  because  of  its  abstractness  it  can  be  used  to  model  many  different 
phenomena.  Every  mathematical  structure  has  an  associated  language  for  mak- 
ing assertions.  If  the  mathematical  model  is  successful,  the  language  of  its 
mathematical  structure  can  be  used  to  make  assertions  about  the  object  being 
modeled. 

The  third  component  of  a  model  is  a  specification  of  the  way  in  which  the  real 
world  is  represented  by  the  mathematical  structure,  that  is,  a  correspondence 
between  the  elements  of  the  first  component  and  those  of  the  second. 

3.  The  Euclidean  Model 

When  using  mathematics  to  describe  man-made  objects,  the  Euclidean 
model  (standard  Euclidean  geometry)  is  usually  satisfactory.  Its  structure  is 
simple  and  pure,  which  appeals  to  an  engineer's  nature.  But  as  technology 
expands  and  we  need  to  describe  processes  that  are  not  well  behaved,  we  need  to 
develop  a  geometry  that  can  adequately  model  our  process  within  a  certain 
closeness  of  scale. 

No  model  can  completely  describe  a  natural  object  because  nature  does 
not  follow  the  man-made  rules  that  we  impose  on  our  model.  But  at  a  given 
scale,  the  model  (if  it  is  accurate)  can  describe  the  object  with  enough  precision 
to  be  of  help  in  constructing  it.  Engineers  use  the  geometry  of  a  straight  line  to 
describe  a  wall  but  this  wall,  when  viewed  closely  enough,  is  not  straight  at  all. 
This  is  of  no  matter  to  the  engineer  because  his  model  is  accurate  for  his  scale  of 
reference. 
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B.  FRACTAL  GEOMETRY 

One  man  who  saw  a  need  for  a  new  geometry  was  Benoit  B.  Mandlebrot.  He 
felt  that  Euclidean  geometry  was  not  satisfactory  as  a  model  for  natural  objects. 
To  anyone  who  has  tried  to  draw  a  picture  of  a  nonregular  object  (such  as  a  tree) 
on  a  computer  graphics  screen,  using  the  Euclidean  drawing  primitives  usually 
provided,  this  is  an  obvious  statement.  The  strength  of  Mandlebrot's  finding  was 
his  research  into  the  findings  of  the  earlier  mathematicians  and  the  development 
of  a  practical  application  of  their  theory.  Mandlebrot  coined  the  term  Fractal  to 
describe  a  class  of  functions  first  discovered  by  Cantor  (Cantor's  dustj.  Koch  (the 
Koch  curve)  and  Peano.  He  showed  how  these  functions  yield  valuable  insight 
into  the  creation  of  models  for  natural  objects  such  as  coastlines  and  mountains. 
Mandlebrot  popularized  the  notion  of  a  fractal  geometry  for  these  types  of 
objects.  Although  he  did  not  invent  the  ideas  he  presents,  Mandlebrot  must  be 
considered  important  because  of  his  synthesis  of  the  theory  at  a  time  when 
science  was  reaching  out  for  new  more  accurate  models  to  describe  its  processes. 

C.  GOALS  OF  THIS  RESEARCH 

There  are  two  approaches  that  can  be  taken  in  the  investigation  of  fractal 
geometry  and  computer  graphics. 

-  To   view    the   computer   as    a   tool    to    enhance    the    investigation    of  fractal 
geometry. 

or 

-  To   view   fractal   geometry   as  a   tool   to   enhance   the   realism  of  computer 
graphics. 

This  research  will  take  the  later  approach1.  It  is  designed  to  investigate  the 
mathematics  of  fractal  geometry  and  to  show  its  application  to  computer 
graphics.  I  hope  to  be  able  to  tame  the  subject  of  fractal  geometry  by  making  its 
mathematics  and  technique  accessible  to  the  average  computer  scientist. 


1  Where  Mandelbrot  took  the  former. 
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II.    THE  MATHEMATICAL  BASIS  OF  FRACTAL  GEOMETRY 

This  chapter  is  a  brief  introduction  to  the  mathematical  foundations  that 
underlie  the  theory  of  fractals.  Little  technique  currently  exists  for  the  practical 
application  to  attain  complete  mathematical  rigor  when  using  fractal  functions 
(i.e.  it  is  very  difficult  to  prove  that  a  set  is  fractal).  This  causes  the  non- 
mathematician  to  accept  much  of  what  he  does  with  fractals  on  faith.  It  is 
instead  important  to  understand  the  theory  intuitively.  This  can  be  gained  by  a 
cursory  look  at  the  mathematical  foundations  for  fractals. 

A.    PRELIMINARIES 

A  complete  definition  of  fractals  is  given  later  in  this  chapter  but  before  we 
can  understand  that  definition,  we  must  establish  a  foundation  in  set  theory. 
Fractals  were  discovered  in  set  theory  and  topology.  They  can  be  considered  as 
an  outgrowth  of  investigations  into  these  related  fields. 

1.  What  is  a  Set 

A  set  is  defined  in  [Ref.  3:pp.  11]  in  the  following  fashion: 

A  set  is  formed  by  the  grouping  together  of  single  objects  into  a  whole.  A  set  is 
a  plurality  thought  of  as  a  unit.  We  can  consider  these  statements  as  expository, 
as  references  to  a  primitive  concept,  familiar  to  us  all,  whose  resolution  into 
more  fundamental  concepts  would  perhaps  be  neither  competent  nor  necessary. 
We  will  content  ourselves  with  this  construction  and  will  assume  as  an  axiom 
that  an  object  M  inherently  determines  certain  other  objects  a,  b,  c,  ...  in  some 
undefined  way,  and  vice  versa.  We  express  this  relation  by  the  words:  The  set  M 
consists  of  the  objects  a,  b,  c,  ... 

This    definition    is    intentionally    vague   to    allow    the    set    to    become    the    basic 
building  block  for  all  mathematical  constructs. 

2.  Some  Set  Theoretic  Concepts 

This  section  presents  some  background  definitions  for  concepts  used  in 
the  body  of  this  chapter.  The  reader  is  directed  to  the  references  for  a  detailed 
explanation  or  proof  [Ref.  2:Chap  5]. 

Definitions: 

Cardinality 

Two  sets  S  and  T  are  said  to  have  the  same  number  of  elements,  or  to  have  the 
same  cardinality,  if  there  is  a  one-one  function  /from  5  to  T. 
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Finite  and  Infinite  Sets 

A  Set  5  is  said  to  be  finite  if  5  has  the  same  cardinality  as  0,  or  if  there  is  a  po- 
sitive integer  n  such  that  S  has  the  same  cardinality  as  {1,2,3,4, 5, ...  n  }.  Other- 
wise S  is  said  to  be  infinite. 

Countability 

A  set  S  is  said  to  be  countable  if  S  has  the  same  cardinality  as  a  subset  of  N,  the 
set  of  positive  integers.    Otherwise,  5  is  said  to  be  uncountable. 

Propositions: 

a).  Any  subset  of  a  finite  set  is  finite. 

b).  Any  subset  of  any  countable  set  S  is  countable. 

c).  The  set  of  natural  numbers  N  is  countable. 

d).  The  set  of  rational  numbers  Q  is  countable. 

d).  The  set  of  irrational  numbers  is  countable. 

e).  The  union  of  a  countable  collection  of  countable  sets  is  countable. 

f).  The  set  of  real  numbers  R  is  uncountable. 

3.     Some  Topological  Concepts 

This  section  presents  some  topological  background  concepts  used  in  the 
body  of  this  chapter.  The  reader  is  directed  to  the  references  for  a  detailed 
explanation  [Ref.  3]. 

Concepts: 

Metric  Space 

A  metric  space  is  a  set  in  which  we  have  a  measure  of  the  closeness  or  proximity 
of  two  elements  of  the  set,  that  is,  we  have  a  distance  defined  on  the  set.  For  ex- 
ample, a  metric  on  R2  would  be  the  pythagorean  metric: 


D((xi,yi),(x2,y2))  =  x/(x1-i2)2+(yi-J/2)2 

Covering  a  Set 

Let  A'  be  a  topological  space  and  S  a  subset  of  X.  A  cover  of  the  set  5  is  exactly 
what  its  name  implies,  a  collection  of  subsets  of  X  which  cover  S,  that  is,  whose 
union  contains  S. 

4.      What  is  a  Function 

A  function  is  defined  in  [Ref.  2:pp.  193-194]  as 

A  function  from  a  set  A  to  a  set  B  is  a  rule  which  specifies  an  element  of  B  for 
each  element  of  A. 

Definition: 

Let  A  and  B  be  sets.  A  function  (or  map,  or  transformation)  f  from  A  to  B. 
denoted  /  :  A  -»  B  .  is  a  relation  from  A  to  B  such  that  for  every  a  which  is  an 
element  of  A,  there  exists  a  unique  b  in  B  such  that  <  a,b>  is  an  element  of  /. 
We  write  f(a)  =  b. 
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If  /is  function  from  A  to  B,  then  A  is  called  the  domain  of  the  function  /and  B 
is  called  the  codomain  of  /. 

To  completely  define  a  function  we  must  specify  the  domain,  the  codomain,  and 
the  value  f(x)  for  each  possible  argument  x. 

Functions  can  be  viewed  as  a  specification  of  a  method  to  describe  the 
creation  of  a  set  from  other  sets  using  some  agreed  upon  mathematical 
symbolism.  The  functions  can  yield  powerful  results  when  the  target  set  (co- 
domain)  is  complex  and  not  easily  described  by  set  theoretical  constructs.  This  is 
especially  true  in  fractal  functions  when  the  domain  is  R1  and  the  object  created 
(set,  co-domain)  is  a  nonregular  shape.  This  is  one  reason  why  the  computer 
graphics  system  is  useful  in  the  investigation  of  fractal  functions.  The  computer 
can  model  the  infinite  function  and  display  a  finite  approximation  of  the  created 
fractal  set. 

5.     Useful  Functional  Concepts 

It  is  often  helpful  to  clearly  understand  the  universe  of  discourse  within 
which  a  function  exists.  The  function  can  be  rigorously  defined  within  the  above 
constructs  but  lack  intuitive  appeal  due  to  its  complexity.  Mathematicians  have 
defined  many  useful  concepts  to  describe  functions.  The  concepts  applicable  to 
this  study  are  described  below. 

a.  Partial  Functions 

Most  of  the  fractal  functions  in  this  study  have  as  their  domain  some 
undefined  subset  of  R  .  It  is  useful  then  to  consider  them  as  partial  functions 
and  not  concern  ourselves  with  a  rigorous  description  of  the  domain  of  the 
fractal.    We  take  our  definition  from  [Ref.  2:pp.  201-202]. 

It  is  often  convenient  to  consider  a  function  from  a  subset  A '  of  A  to  a  set  B 
without  exactly  specifying  the  domain  A  '  of  the  function.  Alternatively,  we  can 
view  such  a  situation  as  one  where  a  function  has  a  domain  A  and  codomain  B, 
but  the  value  of  the  function  does  not  exist  for  some  arguments  of  A.  This  is 
called  a  partial  function. 

Definition: 

Let  A  and  B  be  sets.  A  partial  function  f  with  domain  A  and  codomain  B  is  any 
function  from  A'  to  B,  wnere  A'  is  a  subset  of  A.  For  any  x  which  is  an  element 
of  A  -  A',  the  value  of  /  (x )  is  said  to  be  undef  ined . 

b.  Bijectivity 

It  is  often  useful  to  know  to  what  extent  a  function  maps  from  the 
domain  to  the  codomain.  If  a  function  is  not  a  partial  function  and  every  point  of 
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the  domain  maps  to  a  point  in  the  codomain  then  we  want  to  know  if  all  points 
of  the  domain  A  in  the  mapping  /(A)  map  to  distinct  points  in  the  codomain  B. 
We  may  also  want  to  know  to  what  extent  the  mapping  /(A)  covers  the  set  B. 
The  definition  of  bijective,  surjective  and  injective  functions  is  from  [Ref  2:pp. 
204]. 

Definition: 

Let  /    be  a  function  /:  A  »B. 

(a)  /    is  surjective  (onto)  if  /(A)  =  B, 

(b)  /    is  injective  (one-to-one)  if  a  /  a  '  implies  /(a)  #  /(a'), 

(c)  /    is  bijective  (one-to-one  and  onto)  if  /  is  both  surjective  and  injective. 

6.  Functions  From  RN  -  RN 

A     point     in     R'      space     is     specified     by     an     n- tuple     of     the     form 

(jr1,x2,x3, xn  ).    To  completely  specify  a  function  from  R^  —  RN  each  point  in 

the  domain  must  map  to  a  point  in  the  codomain.    An  example  is: 

/    :  R2-R2 

/((x1?x2))=  (x,2,,/) 

This  function  is  well  defined.  For  each  point  in  the  domain  of  the  function  we 
have  specified  a  unique  point  in  the  codomain. 

Most  of  the  functions  that  are  covered  in  this  study  are  mappings  within 
R2  or  R3.  Fractal  sets  exist  in  all  finite  dimensions  but  it  is  impractical  at  this 
point  to  use  fractal  functions  beyond  the  fourth  dimension  in  view  of  the  graphics 
display  medium's  limitation  to  two  dimensions.  The  use  of  fractal  functions 
whose  dimension  is  between  3  and  4  is  currently  being  investigated  by  allowing 
the  function  to  roam  the  fourth  dimension  and  then  taking  time  slices  which 
yield  three  dimension  approximations  of  the  set  [Ref.  4]. 

7.  Inductive  Definitions  of  Sets 

Functional  constructs  do  not  always  provide  a  convenient  means  of 
charactering  an  infinite  set.  It  is  sometimes  more  eloquent  and  powerful  to  use 
the  inductive  method  to  characterize  a  set. 
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Our  definition  of  inductive  definitions  of  sets  is  from  [Ref.  2:pp.  199-201]. 

An  inductive    definition  of  a  set  always  consists  of  three  distinct  components. 

1.  The  basis,  or  basic  clause,  of  the  definition  establishes  that  certain  objects 
are  in  the  set.  The  basic  clause  establishes  that  a  set  is  not  empty  and  charac- 
terizes the  "building  blocks"  (the  seeds  of  the  induction)  which  are  used  to  con- 
struct the  set  from  the  inductive  clause. 

2.  The  induction,  or  inductive  clause,  of  an  inductive  definition  establishes 
the  ways  in  which  elements  of  the  set  can  be  combined  to  obtain  new  elements. 
The  inductive  clause  always  asserts  that  if  objects  z  ,y  ,...,z  are  elements  of  the 
set.  then  they  can  be  combined  in  certain  specified  ways  to  create  other  elements 
of  the  set  (thus  from  the  basic  clause  (or  seeds)  of  the  induction  we  induce  the 
remaining  elements  of  the  set). 

3.  The  extremal  clause  asserts  that  unless  an  object  can  be  shown  to  be  a 
member  of  a  set  by  applying  the  basis  and  inductive  clauses  a  finite  number  of 
times,  then  the  object  is  not  a  member  of  the  set. 

An  example  of  an  inductively  defined  set  is: 

(Basis) 

0  t  A 

(Induction) 

If  n   e  A,  then  (n  +2)  e   A 

(Extremal) 

No  integer  is  an  element  of  A  unless  it  can  be  shown  to  be  so  in  a  finite  number 
of  applications  of  clauses  1  and  2  above. 

The  set  that  we  defined  is  the  set  of  all  even  nonnegative  integers. 

8.     The  Path  To  Fractals 

The  path  to  fractals  by  the  non-mathematician  is  not  through  theory  but 
through  the  investigation  of  their  functions  and  methods  of  construction.  This 
investigation  (and  experimentation)  yields  considerable  insight  into  the  nature  of 
fractals. 

The  choice  of  which  set-descriptive  methodology  to  use  (functional  or 
inductive)  in  describing  a  set  is  often  a  matter  of  style  but  can  be  dictated  by 
necessity  if  one  method  is  inordinately  tedious. 

Most  of  the  fractal  functions  that  are  introduced  in  this  study  use  the 
inductive  method  as  the  primary  functional  tool.  In  fact,  these  functions  are  a 
hybrid  of  the  functional  and  inductive  constructs  described  above. 


Often  called  a  recurrence  definition. 
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B.     DIMENSION 

The  classification  of  fractal  sets  from  non-fractal  sets  is  based  on  the 
dimensional  qualities  of  the  set.  To  understand  fractals  you  must  have  an 
appreciation  for  these  differences. 

The  concept  of  dimension  is  one  rife  with  difficulties.  Many  of  the  gr^at 
mathematicians  have  attempted  to  define  dimension  as  a  rigorous  concept 
consistent  with  the  known  mathematical  systems.  For  each  of  their  attempts 
however,  the  concept  becomes  more  prone  to  contradiction  and  paradoxes. 

There  currently  exist  five  definitions  of  dimension  that  date  back  to  the  late 
1800V.  The  classification  of  Fractal  sets  into  a  class  of  sets  is  the  result  of  the 
discovery  of  functions  that  created  sets  which  did  not  fit  comfortably  into  the 
topological  definition  of  dimension  (which  was  the  accepted  definition  at  the 
time).  Fractal  sets  are  rigorously  classified  as  those  sets  that  demonstrate  a 
difference  between  the  Hausdorff-Besicovitch  dimension  and  the  standard 
topological  dimension. 

1.     An  Intuitive  Approach  to  Dimension 

Dimension  is  a  concept  that  seems  intuitive  when  it  is  first  introduced  in 
Euclidean  geometry  as  the  standard  three  dimensions.  Long  after  Euclid  made 
the  first  attempts  at  defining  dimension  and  concurrent  with  the  discovery  of 
atomic  particle  physics  the  concept  of  dimension  was  rethought  by  the  prominent 
mathematicians  of  the  time.  This  was  necessary  to  realign  the  mathematical 
model  for  the  geometry  of  objects  with  the  new  view  of  what  those  objects  were 
made  of.  Our  increasing  ability  to  focus  on  the  nature  of  matter  inevitably 
causes  the  models  we  use  to  change. 

The  dilemma  that  arose  from  the  new  concepts  of  dimension  quickly 
developed  into  a  theoretical  debate  that  left  intuition  behind.  When  human 
intuition  fails,  we  must  rely  upon  well  founded  models  that  are  based  on  axioms 
of  basic  mathematical  truth.  It  is  only  through  the  rigorous  investigation  of  our 
mathematical  models  that  allows  us  to  go  beyond  intuition  and  investigate  the 


1  1).  Cantor  and  Minkowski;  2).  Bouligand  and  Minkowski;  3).  Pontrajgin.  Schnirelman  and 
Kolornogorov,  Tihomirov;  4).  Hausdorff-Besicovitch  and  5).  the  topological  dimension  (there  are 
others).  Most  of  these  definitions  are  concerned  with  the  most  efficient  method  of  covering  a  set 
(i.e.  are  topological  concerns). 
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true  physical  nature  of  the  objects  we  model.    The  debate  still  rages  today  and 
borders  on  the  philosophical.    Two  examples  should  suffice  to  demonstrate  the 
complexity  and  possible  paradoxes  that  can  arise  from  dimension  theory. 
The  first  is  from  [Ref.  5:pp.  323-344]. 

Consider  the  way  in  which  we  define  the  density  of  air  at  a  given  point  and  at  a 
given  moment.  We  picture  a  sphere  of  volume  V  centered  at  that  point  and  in- 
cluding the  mass  Wi.  The  quotient  M/V  is  the  mean  density  within  the  sphere, 
and  by  the  true  density  we  denote  some  limiting  value  of  this  quotient.  This  no- 
tion, however,  implies  that  at  the  given  moment  the  mean  density  is  practically 
constant  for  spheres  below  a  certain  volume.  This  mean  density  may  be  notably 
different  for  spheres  containing  1,000  cubic  meters  and  1  cubic  centimeter 
respectively. 

Suppose  the  volume  becomes  continually  smaller.  Instead  of  becoming  less  and 
less  important,  these  fluctuations  come  to  increase.  For  scales  at  which  the 
brownian  motion  shows  great  activity,  fluctuations  may  attain  1  part  in  1,000, 
and  they  become  of  the  order  of  1  part  in  5  when  the  radius  of  the  hypothetical 
spherule  becomes  of  the  order  of  a  hundredth  of  a  micron. 

One  step  further  and  our  spherule  becomes  of  the  order  of  a  molecule  radius.  In 
a  gas  it  will  generally  lie  in  intermolecular  space,  where  its  mean  density  will 
henceforth  vanish.  At  our  point  the  true  density  will  also  vanish.  But  about 
once  in  a  thousand  times  that  point  will  lie  within  a  molecule,  and  the  mean 
density  will  be  a  thousand  times  higher  than  the  value  we  usually  take  to  be 
the  true  density  of  the  gas. 

Let  our  spherule  grow  steadily  smaller.  Soon,  except  under  exceptional  cir- 
cumstances, it  will  Decome  empty  and  remain  so  henceforth  owing  to  the  intra- 
atomic  space;  the  true  density  vanishes  almost  everywhere,  except  at  an  infinite 
number  of  isolated  points,  where  it  reaches  an  infinite  value. 

The  second  is  from  [Ref.  l:pp.  17-18]. 

Consider  a  ball  of  10  cm  diameter  made  of  a  thick  thread  of  1  mm  diameter  that 
(in  latent  fashion)  possesses  several  distinct  effective  dimensions. 

To  an  observer  placed  far  away,  the  ball  appears  as  a  zero-dimensional  figure:  a 
point.  As  seen  from  a  distance  of  10  cm  resolution,  the  balj  of  thread  is  a  three- 
dimensional  figure.  At  10  mm,  it  is  a  mess  of  one-dimensional  threads.  At  0.1 
mm,  each  thread  becomes  a  column  and  the  whole  becomes  a  three-dimensional 
figure  again.  At  0.01  mm,  each  column  dissolves  into  fibers,  and  the  ball  again 
becomes  one-dimensional,  and  so  on,  with  the  dimension  crossing  over  repeated- 
ly from  one  value  to  another.  When  the  ball  is  represented  by  a  finite  number  of 
atomlike  pinpoints,  it  becomes  zero-dimensional  again. 

It  is  interesting  to  note  that  each  of  these  examples  demonstrate 
dimension  as  a  reflection  of  physical  properties  dependent  on  the  observers  point 
of  reference.  Each  ends  with  reference  to  the  paradox  of  atomic  particles.  That 
paradox  is,  for  any  collection  of  finite  (or  countably  infinite)  points,  the 
dimension  is  zero  [Ref.  6:pp.  1-8].  Since  the  earth  and  sun  each  have  a  finite 
collection  of  atoms  then   accordingly  their  dimension  is  zero.   Dimension  exists 
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only  for  a  mathematical  continuum  and  as  such  lacks  application  to  the  physical 
universe  as  we  currently  know  it4. 

The  two  properties  (continuity  and  dimension)  cannot  be  separated. 
Before  the  advent  of  atomic  theory,  matter  was  viewed  as  continuous  and 
composed  of  basic  elements  that  were  indivisible.  While  the  debate  raged  over  the 
practical  and  philosophical  aspects  of  the  nature  of  matter  it  became  apparent 
that  the  mathematical  models  which  represent  matter  would  have  to  change.  It  is 
not  practical  to  represent  objects  by  representing  each  atom  and  its  position 
relative  to  the  entire  set.  The  power  of  modeling  would  thus  be  lost;  that  is.  the 
ability  to  model  complex  objects  and  their  interaction  by  relatively  simple 
constructs.  Thus  the  fact  is  reinforced  that  models  can  only  represent  objects 
through  gross  approximations  and  that  the  model  is  only  effective  for  a  restricted 
frame  of  reference.  Without  this  realization,  dimension  would  have  very  little 
application. 

2.      Topological  Dimension 

a.     An  Intuitive  Approach 

The  concept  of  dimension  is  very  old.  It  is  based  on  the  algebraic 
concepts  of  Euclidean  n  space  and  the  notion  that  a  set  has  dimension  n  if  the 
least  number  of  real  parameters  needed  to  describe  its  points  was  n.  This  fuzzy 
definition  was  accepted  for  a  very  long  time  until  the  advent  of  Cantor  and  set 
theory.  Cantor  showed  that  dimension  can  be  changed  by  a  1-1  transformation 
from  an  interval  to  a  planar  object.  The  fuzzy  notion  of  dimension,  as  defined, 
was  challenged  and  required  rethinking. 

The  mathematicians  who  did  not  accept  many  of  the  findings  of  set 
theory  at  the  time  (but  who  could  not  disregard  Cantors  findings)  began  to 
consider  ways  of  explicitly  defining  dimension.  The  new  definition  would  have  to 
be  applicable  to  the  bizarre  functions  of  Cantor.  Koch  and  Peano  as  well  as  the 


4  The  set  theoretical  concepts  of  finiteness,  countably  infinite  and  uncountabl)  infinite  (con- 
tinuous) sets  carry  with  them  very  profound  implications.  It  is  premature  to  view  matter  a>  mere- 
ly collections  of  finite  atoms.  Science  may  yet  find  true  continuity  (in  the  mathematical  sense)  in 
atomic  matter  and  the  universe.  For  now,  matter  is  what  it  is  and  our  pronunciations  upon  it  will 
not  change  its  true  texture. 
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relatively  simple  objects  that  had  previously  fit  into  the  old  definition  without 
contradiction. 

There  was  one  crucial  problem  that  Cantor's  findings  raised 
[Ref.  6:pp.  4-6]: 

An  extremely  important  question  was  left  open:  Is  it  possible  to  establish  a 
correspondence  between  Euclidean  n  space  ana  Euclidean  m  space  combining  the 
features  of  both  Cantor's  and  Peano's  constructions,  i.e.  a  correspondence  which 
is  both  1:1  and  continuous?  The  question  is  crucial  since  the  existence  of  a 
transformation  of  the  stated  type  between  Euclidean  n  space  and  Euclidean  m 
space  would  signify  that  dimension  (in  the  natural  sense  that  Euclidean  n  space 
has  dimension  n  )  has  no  topological  meaning  whatsoever. 

This  fundamental  problem  was  answered  in  1911  by  Brouwer.  He 
proved  that  Euclidean  n  space  and  Euclidean  m  space  were  not  homeomorphic 
unless  n  equals  m.  To  say  that  two  spaces  A  and  B  are  homeomorphic  means 
that  a  mapping  /  :A^B  exists,  such  that  /is  continuous  over  .4  and  bijective. 
Additionally,  the  inverse  of  this  mapping  /  ~l  B^A  ,  is  continuous  over  B  and 
bijective.  If  two  spaces  are  homeomorphic  then  it  is  analogous  to  saying  that  they 
are  topologically  equivalent. 

Further   research   was   done   and   a  precise  definition   of  topological 
dimension  of  a  set  was  developed.  This  definition  assigned  an  integer  value  as  the 
dimension  of  any  set  based  on  its  topological  properties, 
b.     Definition  of  Topological  Dimension  n 

The  rigorous  investigation  of  dimension  is  beyond  the  scope  of  this 
study  but  the  following  definition  is  included  for  completeness  [Ref.  6:pp.  24]: 

Roughly  speaking,  we  may  say  that  a  space  has  dimension  ^  n  if  an  arbitrarily 
small  piece  of  the  space  surrounding  each  point  may  be  delimited  by  subsets  of 
dimension  <  n  -  1.  This  method  of  definition  is  inductive,  and  an  elegant  start- 
ing point  for  the  induction  is  given  by  prescribing  the  null  set  as  the  (-1)  dimen- 
sional space. 

Definition: 

The  empty  set  and  only  the  empty  set  has  dimension  -1. 

A  space  A'  has  dimension  ^  n  (  n  >  0)  at  a  point  p  if  p  has  arbitrarily  small 
neighborhoods  whose  boundaries  have  dimension  ^  n  -  1. 

X  has  dimension  ^  n,  dim  X  ^  n,  if  .V  has  dimension  <  n  at  each  of  its  points. 

X  has  dimension  n  at  point  p  if  it  is  true  that  X  has  dimension  -^  n  at  p  and  it  is 
false  that  X  has  dimension  ^  n  -  1  at  p. 

X  has  dimension  n  if  dim  X  <  n  is  true  and  dim  X  ^  n  -  1  is  false. 

X  has  dimension  oo  if  dim  X  ^  n  is  false  for  each  n. 
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The  topological  dimension  is  rigorous  and  consistent  for  all  sets  that 
exist  within  a  metric  space.  The  problem  that  arises  with  fractal  sets  and  its 
topological  dimension  is  not  that  the  topological  dimension  is  wrong.  Fractal  sets 
like  all  sets  in  a  metric  space  exhibit  a  topological  dimension.  The  question  is 
then,  is  the  topological  dimension  an  accurate  description  of  the  dimension  of  the 
set  or  can  we  find  a  better  way  to  characterize  the  dimensional  qualities  of  the 
set?  This  question  can  be  extended;  is  the  topological  definition  of  dimension 
useful  and  consistent  with  the  notion  of  dimension  and  space?  Can  we  devise  a 
better  concept  which  can  further  refine  dimension  and  make  it  more  useful3. 
3.     The  Hausdorff-Besicovitch  Dimension 

This  section  is  intentionally  brief  due  to  the  subject's  complexity  and  to 
the  lack  of  practical  technique  that  it  yields.  The  Hausdorff  measure  of  a  set  is  a 
complex  characterization  of  a  method  for  covering  a  set.  Hausdorffs  theorem  is 
proved  using  the  existential  qualities  of  infinite  sets  in  a  metric  space  H.  While 
the  theorem  may  be  important  to  mathematical  theory,  it  proves  unfortunate 
that  there  is  no  straightforward  practical  method  for  determining  the  Hausdorff 
measure  of  a  set. 

a.     An  Intuitive  Approach  to  the  Hausdorff  Dimension 

The  acceptance  of  the  Hausdorff  method  for  covering  a  set  as  a 
measure  of  dimension  is  not  universal  [Ref.  6:pp.  102]  and  [Ref.  l:pp.  363-365]. 
The  debate  is  between  the  disciplines  of  topology  and  metrics  and  is  not  wholly 
germane  to  this  study.  It  is  beneficial  to  divorce  ourselves  from  the  debate  and 
consider  both  the  topological  dimension  DT  and  the  Hausdorff-Besicovitch 
dimension  (HB)  as  merely  measures  of  qualities  of  a  set's  structure.  Certainly 
sets  exist  that  have  a  topological  dimension  equal  to  1  but  in  no  way  resemble  a 
simple  Euclidean  curve  If  the  Hausdorff  dimension  yields  a  better  measure  of  a 
set's  structure  that  provides  a  mathematical  and  intuitive  difference  that  is  useful 
to  us,  it  would  be  beneficial  for  us  to  investigate  it. 


5  Try  not  to  ascribe  grandiose  implications  to  a  set's  dimension  (the  fourth  dimension  as  time 
or  some  such)  as  this  is  premature  at  best.  Rather,  view  a  set's  dimension  as  merely  descriptive 
terminology  much  like  the  terminology  of  bijectivity  describes  a  function's  characteristics  The 
problem  most  people  have  with  this  mental  abstraction  is  the  visual  reinforcement  that  the\  re- 
ceived from  the  notion  of  the  standard  three  dimensions. 
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The  Hausdorff  measure  of  a  set  was  developed  during  the  same 
period  that  the  new  topological  dimension  was  invented  to  solve  the  paradoxes  of 
Cantor.  The  topological  dimension  was  based  on  the  idea  of  a  neighborhood  of  a 
point  within  a  Euclidean  space  of  R  .  The  connection  to  metric  spaces  and  the 
idea  of  measure  is  obvious  when  you  consider  that  the  Hausdorff  measure  of  a  set 
is  also  based  on  this  notion  of  a  spherical  neighborhood  and  what  Hausdorff  calls 
the  test  function  of  a  set.  The  test  function  of  a  set  denoted  h(p  )  is  a  function 
that  characterizes  the  "best"  method  of  covering  a  point  with  the  spherical  ball 
of  radius  p  that  covers  points  of  the  set,  which  in  their  union,  cover  the  entire 
set. 

Consider  for  example  the  test  function  for  a  surface  within  R\  A 
surface  can  be  covered  by  discs  (circles).  The  formula  for  the  area  of  a  circle 
becomes  the  test  function  for  the  surface.  The  formula  for  the  area  of  a  circle 
always  contains  the  constant  factor  it  multiplied  by  the  square  of  the  radius  r. 
This  radius  is  the  measure  p  as  above.  This  leaves  us  with  a  test  function  for  a 
planar  shape  in  R2  of  h(p  )  =  wp2. 

You  might  expect  that  the  test  function  for  a  spherical  neighborhood 
in  a  Euclidean  space  above  R3  would  be  very  difficult  to  imagine,  and  indeed  it 
is.  Hausdorff  further  complicated  the  idea  of  test  functions  (even  within  the  lower 
dimensions)  by  allowing  a  test  function  to  assume  a  non-integer  parameter  d  so 
that  the  test  function  h(p  )  =  0(d)p  could  have  a  real-valued  parameter  d. 
This  further  refinement  of  the  test  function  allowed  Hausdorff  to  make  assertions 
about  how  this  test  function  h(p  )  behaved  when  the  parameters  p  and  d  were 
allowed  to  vary. 

Hausdorff  imagined  the  parameter  p  reducing  in  size  until  it 
approached  zero.  The  effect  of  this  on  our  disc  example  is  increasingly  smaller 
and  smaller  discs  around  points  of  the  planar  set.  As  the  disc  size  is  decreased, 
fewer  points  of  the  set  are  contained  in  each  disc  neighborhood.  This  requires 
more  discs  to  cover  the  set.  As  the  parameter  p  becomes  infinitely  small  the 
number  of  discs  required  to  cover  the  planar  set  approaches  oc.  We  allow  the 
parameter  to  grow  arbitrarily  small.    It  is  interesting  to  study  the  test  function 
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and  see  what  happens  to  the  total  area  when  the  areas  of  the  collection  of  discs 
which  cover  the  planar  set  are  summed. 

Let's  reflect  upon  the  mathematical  process  that  we  are  developing. 
When  we  attempt  to  approximate  the  area  of  the  planar  set  by  the  union  of  the 
discs  which  make  up  its  cover,  we  are  essentially  observing  small  patches  of  the 
surface  and  approximating  the  area  of  the  set  by  making  assertions  about  the 
intrinsic  qualities  of  the  patch.  The  notion  that  this  measure  is  merely  an 
approximation  is  important.  As  the  size  of  our  patch  grows  increasingly  small,  we 
can  expect  that  we  will  get  a  better  fit  with  our  patches  and  hence  a  better 
approximation  of  the  area.  The  notions  of  approximation  and  fit  become 
increasingly  helpful  when  you  consider  functions  which  describe  sets  of  infinitely 
rough  texture  as  we  find  in  fractal  sets. 

The  importance  of  Hausdorffs  discovery  lies  in  the  fact  that  for  a 
test  function  of  a  set,   the  parameter  d  is  special.     As  p-»0  he  discovered  that 

there  existed  a  unique  real  number  d  such  that  for  d'<d.  the  infinum  defined 
by  the  test  function  using  d'  approaches  oc  (for  any  countably  infinite  set).  And 
for  a  d'>d  the  corresponding  infinum  approaches  zero.  This  means  that  every 
set  has  a  parameter  which  can  be  associated  with  it  that  is  closely  related  to  the 
amount  of  space  that  it  occupies.  This  number  is  the  Hausdorff-Besicovitch 
dimension  (measure)  of  the  set. 

These  results  only  tell  us  that  a  number  and  function  exist.  They  do 
not  tell  us  how  to  compute  them  in  the  general  case.  This  gives  us  the  quandary 
of   dealing    with    the    HB    dimension    as    a    known    concept    that    we    can    make 
allusions  to,  but  can  rarely  compute  (at  least  at  the  present  time), 
b.     Definition  of  the  Hausdorff  Dimension 

The  formal  definition  of  the  Hausdorff  dimension  requires  that  we 
formalize  the  inttiitive  discussion  above.  We  first  define  what  a  ;>-measure  of  a  set 
is  (analogous  to  the  disc  above)  [Ref.  6:pp.  102-103]. 

A  p-dimensional  measure  for  each  non-negative  real  number  p  was  defined  bv 
Hausdorff  for  arbitrary  metric  spaces.  This  measure  is  a  metrical  concept.  whiTe 
dimension  is  purely  "topological.  Nevertheless  there  is  a  strong  connection 
between  the  two  concepts,  for  it  turns  out  that  a  space  of  [topological]  dimen- 
sion n  must  have  positive  rt-dimensional  measure  (the  Hausdorff  measure). 


For  our  example  this  would  be  the  sum  of  the  areas  of  the  discs 
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Definition: 

Let  X  be  a  space  and  p  an  arbitrary  real  number.  0  ^  p  <oo.  Given  e  >0  let 

m<p   =  inffj  [6{A,)\> 


1  =  1 


where  X  =  A1  +  A2  +  A*  +  •  •  ■  is  any  decomposition  of  X  into  a  countable 
number  of  subsets  of  diameter  less  than  e  ,  and  the  superscript  p  denotes  ex- 
ponentiation. Let 

m.  (X)  =  sup  m%  (X). 
t  >o 

mp  (X)  is  called  the  p-measure  or  (p-dimensional)  measure  of  X. 

Proposition: 

If  p  <q  then  mp  (X)  ^  m?  (X);  in  fact  p  <q  and  mp  (X)  <oo  imply  m?  (X)  =  0. 

Conversely,  if  p  >q  then  mp  (X)  ^  mq  (X);  and  if  p  >g  and 
mp  (X)  <oc   implies  mq  (X)  =  oo  . 

c.      Mandelbrot's  Misgivings  about  the  HB  Measure 

It  is  clear  from  the  previous  definition  that  for  practical  applications 
the  Hausdorff  dimension  is  difficult  to  compute  directly.  In  [Ref.  l:pp.  14-19]  and 
[Ref  7],  Mandelbrot  expressed  a  distaste  for  the  focusing  of  attention  on  the  HB 
dimension.  He  states: 

''I  developed  the  definition  of  fractals  using  the  topological  and  Hausdorff  di- 
mensions in  response  to  colleagues  who  urged  me  to  do  so.  They  felt  that  it  was 
necessary  to  rigorously  define  the  concept  within  firm  mathematical  criteria.  I 
have  come  to  believe  that  an  empirical  definition  would  be  more  beneficial  at 
this  time  because  the  present  definition  denies  the  inclusion  of  some  shapes  that 
could  best  be  described  as  fractals." 

Mandelbrot  believes  that  the  definition  of  the  Hausdorff  dimension  is 
too  difficult  to  deal  with  and  is  perhaps  too  restrictive.  He  prefers  to  focus  on  the 
behavioral  aspects  of  the  recurrence  relationship  involved  in  fractals  and  the 
empirical  results  from  these  equations. 

A  practical  application  of  fractal  functions  in  computer  graphics 
does,  by  necessity,  bend  to  this  same  paradigm.  This  realization  should  not  blind 
us  to  the  fundamental  nature  of  a  fractal  equation's  uniqueness,  however.  It  is 
important  to  understand  the  dimensional  aspect  of  the  fractal  discourse  to 
appreciate  the  potential  importance  of  fractals. 

It  is  not  preordained  that  fractal  equations  model  nature  with  a 
greater  degree  of  accuracy  then  does  the  Euclidean  model.  The  future  may  prove 
the  fractal  model  the  superior  method,  however.  The  dimensional  qualities  of 
fractal  functions  may  be  the  aspect  that  proves  this  to  be  so. 
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4.     Why  Consider  Dimension 

For  the  purpose  of  this  study,  it  is  not  important  that  a  rigorous  feel  for 
the  mathematical  properties  of  dimension  theory  be  grasped.  In  fact  one  needs  no 
knowledge  of  dimension  to  use  fractal  techniques  in  the  generation  of  computer 
graphics  terrain.  The  literature  is  rife  with  articles  about  fractal  objects  in 
computer  graphics  and  it  seems  de  rigueur  to  include  an  approximate  fractal 
dimension  as  part  of  its  description.  The  techniques  used  to  approximate 
dimension  as  presented  in  [Ref.  l:pp.  56-57]  are  mathematically  unproven'.  More 
importantly,  the  dimension  yields  little  intuitive  insight:  one  is  hard  pressed  to 
describe  the  differences  between  an  object  with  an  approximate  dimension  of  2.37 
and  another  with  a  dimension  of  2.45.    The  pictures  are  much  more  descriptive. 

One  use  of  approximate  fractal  dimension  is  to  describe  an  object's 
relative  roughness.  It  is  beneficial  to  view  a  fractal  dimension  as  degrees  of 
roughness  between  the  standard  three  dimensions.  If  the  dimension  is  between  1 
and  2  then  the  object  should  be  a  very  irregular  curve.  If  the  dimension  of  that 
curve  approaches  1  then  the  curve  is  probably  not  very  rough  and  would  lack  any 
interesting  diversion  from  an  ordinary  plane  curve.  If  the  dimension  of  the  curve 
approaches  2  then  the  curve  becomes  like  a  plane  or  filled  polygon  and  again 
lacks  appeal.  The  most  interesting  fractal  curves  are  those  which  demonstrate 
dimensions  nearer  the  center  of  the  scale  between  the  standard  Euclidean 
dimensions.  A  similar  argument  can  be  made  for  solid  objects  with  dimensions 
between  2  and  3. 

This  is  not  to  say  that  fractal  geometry  is  not  a  powerful  tool  for  the 
graphics  programmer.  The  evidence  of  the  power  of  fractals  to  model  objects  of 
considerable  complexity  is  clearly  demonstrated.  To  date,  this  power  has  not 
been  matched  by  other  standard  methods. 

The  graphics  programmer  should  not  concern  himself  greatly  with  the 
dimension  that  is  demonstrated  at  different  levels  of  object  construction.  He  must 
concern  himself  with  the  techniques  of  construction  and  the  realism  that  is 
achieved. 


7  Mandelbrot's  method  for  estimating  the  Hausdorff- Be.sicovitch  dimension  for  non-random 
sets  built  through  self  similar  shapes  will  be  introduced  in  section  C  after  the  hoch  curve  i> 
described. 
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C.    FRACTAL  CURVES  AND  SETS 
1.     Definition  of  Fractal  Sets 

We  take  our  definition  of  fractal  sets  from  [Ref  l:Chap  3  and  pp.  361]. 

Definition: 

A   fractal  set   is   a   set   for   which    the    Hausdorff Besicovitch  dimension   strictly 
exceeds  the  topological  dimension. 


As  we  have  established  and  is  emphasized  by  Mandelbrot,  this  definition 
is  not  very  useful.  The  definitions  of  topological  and  Hausdorff  dimensions  are 
very  involved.  It  is  a  gargantuan  effort  to  prove  that  a  set  has  a  topological  or 
Hausdorff  dimension  (if  one  desires  complete  rigor,  typically  the  topological 
dimension  is  derived  by  the  least  parameter  approach  (section  2)).  When  using 
fractal  functions  then,  it  is  practical  to  assume  that  because  the  functional 
method  you  use  is  fractal-like  that  the  dimension  is  fractal. 

The  functional  techniques  to  be  introduced  have  a  certain  methodology 
that  creates  fractal  sets  with  a  behavior  that  is  disciplined  and  predictable.  The 
assumption  is,  since  these  methodologies  are  well  behaved,  that  any  set  created 
by  these  methods  (with  some  careful  restrictions)  will  itself  be  a  fractal  set. 

WTe  are  left  with  a  practical  methodology  whereby  we  discover  fractal 
functional  methods,  prove  that  the  set  created  is  fractal,  characterize  the  fractal 
part  of  the  functional  method   (carefully)   and  then   enshrine  that  method  as  a 

fractal  method.  If  one's  purpose  is  a  practical  application  of  fractal  techniques8 
and   not    a  rigorous  mathematical   investigation   then   this   is   a   reasonable   and 
practical  approach.    This  approach  is  taken  in  the  remainder  of  this  study. 
2.     Constructing  Fractal  Sets 

In  order  to  describe  the  construction  of  the  Koch  curve,  it  is  necessary  to 
present  terminology  introduced  by  Mandelbrot  [Ref.  l:pp.  34-35]. 

We  use  a  geometrical  shape  (at  first  a  straight  line)  and  call  this  shape 
an  INITIATOR.  We  create  another  shape  that  is  constructed  with  shapes 
similar  to  the  initiator  and  call  this  a  GENERATOR.     We  define  a  sequence  of 


8  A  working  model  or  equation  where  you  are  only  concerned  about  the  behavior  and  not  the 
exact  mathematical  properties. 
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transformations  upon  all  current  initiators  (suppose  there  are  m  such  initiators) 
in  the  construction  by  applying  the  generator  to  all  initiators.  This  creates  a  new 
construction  that  consists  of  m  x  r  (where  r  is  the  number  of  distinct  parts  of  the 
generator)  sides  where  each  side  is  a  shape  that  is  similar  to  the  initiator.  The 
next  step  is  to  again  apply  the  generator  to  all  initiators. 

This  recursive  definition  has  no  terminating  event  but  is  continued  ad 
infinitum.  This  functional  process  is  well  suited  for  recursion  because  the 
application  of  the  generator  to  the  initiator  is  constant  with  respect  to  method 
and  varies  only  to  scale.  It  is  also  well  suited  for  parallel  processing  (in  the 
computer  science  sense)  because  each  application  of  a  generator  on  all  current 
initiators  is  independent. 

These  concepts  are  probably  confusing  at  this  point  and  were  especially 
difficult  to  visualize  when  they  were  first  envisioned  because  the  authors  had  few 
tools  beyond  mental  imagery  to  convey  their  point.  This  is  probably  why  they 
were  largely  ignored  for  70  years.  It  is  much  easier  to  visualize  these  functions 
when  they  are  shown  on  a  computer  graphics  display. 
3.     The  Koch  curve. 

The  mathematicians  Koch  and  Cantor  developed  functions  which 
attempted  to  challenge  the  mathematical  models  of  continuity  and 
differentiability.  These  equations  were  developed  during  the  great  debate  on  set 
theory  and  were  used  by  Cantor  in  arguing  for  his  theory.  These  functions  were 
like  none  before,  using  constructions  which  played  upon  natural  geometric 
constructs  but  when  combined  with  the  power  of  infinite  recursion  became  sets 
which  defied  intuition.  It  was  not  until  much  later  that  mathematicians  were  able 
to  reconcile  these  functions  with  algebra  and  set  theory.  The  function  to  be 
introduced   has   been   proven   to   have   a   Hausdorff-Besicoritch   dimension  which 

exceeds  its  topological  dimension  . 

The  Koch  curve  is  a  very  beautiful  curve  that  at  first  gives  the  observer 
the  impression  of  a  snowflake  or  a  coastline  (Fig  2.1).  Mandelbrot  uses  this  type 
of  construction  (with  variation  —  i.e.  randomness  or  less  behaved  generators)  to 
draw  coastlines   that   look  very   realistic.   To  understand   this  construction   is   to 


Thus  fractals  do  exist  and  it  is  possible  to  rigorously  prove  it  to  be  so. 
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understand  one  method  of  obtaining  a  fractal  set  from  a  well  defined  non-fractal 
set  (R2)  using  non-random  techniques.  This  method  of  construction  (in  the 
general  sense)  is  very  powerful  and  is  used  throughout  this  thesis. 

To  construct  the  Koch  curve,  we  use  three  initiators  (line  segments)  of 
equal  length  and  join  them  to  form  an  equilateral  triangle10.  To  construct  the 
generator  we  use  four  line  segments  that  are  each  —  the  length  of  the  initiators 

and  apply  these  to  each  initiator  (Fig  2.1).  This  yields  a  new  geometric  figure 
with  12  sides  versus  the  original  3  and  a  total  perimeter  length  of  4  units  of 
length  versus  the  original  3  units. 

Figure  2.1  demonstrates  the  first  and  second  recursive  iterations  of 
building  the  Koch  curve.  Observe  how  the  progression  develops  to  yield  the  final 
figure  in  Figure  2.1.    Imagine  this  progression  occurring  indefinitely. 

-  With  each  iteration  of  applying  the  generators  the  total  perimeter  length 
increases  by  —  over  the  previous  perimeter. 

-  The  length  of  the  curve  begins  to  increase  without  bound  even  though  the 
length  of  the  initiator  decreases  to  an  infinitely  small  length.  Hence  the 
curve's  length  is  unbounded  with  no  point  intersecting  but  yet  is  contained 
in  a  small  bounded  two  dimensional  area. 

-  The  points  of  the  curve  are  by  construction  only  the  end-points  of  each 
initiator  and  each  point  is  clearly  distinct  from  the  other  (  no  two  points  are 
connected). 

-  Although  each  point  is  distinct  at  any  one  level  of  the  curve  construction,  it 
can  be  proven  that  the  curve  when  viewed  in  the  limit  is  continuous  at  every 
point. 

-  That  due  to  the  above  qualities  the  curve  is  not  differentiable  at  ANY  point. 

It  is  important  to  realize  that  the  endpoints  of  the  lines  (initiators)  are  the  only 
points  of  the  curve.  The  line  only  serves  as  a  vehicle  by  which  the  points  may  be 
easily  determined.  The  exact  same  set  could  be  built  using  180°  arcs  as  initiators. 
An  algorithm  and  computer  graphics  program  for  the  construction  of  the 
Koch  curve  is  presented  in  Chapter  4  and  Appendix  A. 


The  choice  of  an  equilateral   triangle  was  arbitrary.  We  could  have  chosen  any  shape  as 
long  as  it  was  made  up  of  Initiators  and  avoided  intersecting  lines  during  recursion. 
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r  =  - 


THE  KOCH  CURVE 


Recursive  term- 
ination distance 
is  .05  inch. 


DT  =  1 , HB  =  1.26 


Output  medium: 
Laser  printer 

Resolution : 
300  dots/inch 


Figure  2.1  The  Koch  Curve 


27 


4.     Mandelbrot's  Dimension  Approximation  Function 

In  the  above  section,  one  functional  construction  technique  for  building  a 
fractal  set  has  been  introduced.  It  is  possible  to  approximate  the  dimension  of  a 
fractal  curve  that  is  built  using  these  constructs.  In  [Ref.  l:pp  56-57],  Mandelbrot 
introduces  a  function  that  is  based  on  the  similarity  properties  of  the  above 
technique.      This    function    has    a    real   exponent    D    that    is    the    approximate 

dimension  of  the  fractal  set11. 

Consider  a  method  of  paving  (covering?)  a  Euclidean  shape.  Divide  a  line 
segment  into  N  segments  with  each  segment  a  part  of  the  original  segment  such 
that  the  sum  of  the  lengths  of  the  N  segments  equals  the  length  of  the  original 
segment.   It   follows   then,   that   the  sum  of  the  ratios  of  the   divided  segments 


lengths    to    the    original    segment    length 
equal  1. 


/  .   .  \ 

Divided  segment   s 

e.g.   rs  =  — — ; 

Original   segment 


must 


N 

£rs=  1 
s=l 

We  know  that  the  dimension  of  a  line  is  equal  to  1.    If  we  raise  each  of  the  above 

ratios  to  the  power  D  (where  D  =  1)  the  equivalence  still  holds. 

Eh     =1 

Let's  allow  the  Koch  function  to  assume  a  similar  dimensional  relationship  but 
treat  D  as  a  real  valued  unknown.  Refer  once  again  to  Figure  2.1.  Notice  that 
the  length  of  the  four  line  segments  that  make  up  the  generator  have  a  length 

ratio  of  —  to  the  initiator.  Call  this  ratio  rm.  Notice  that  the  generator  is  made 

up  of  four  initiator  shapes.  Call  this  number  N. 


s  M 


=    1 


11    CAUTION:    This    technique    does    not    necessarily    apply    to    other    fractal    functional 
methods. 
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Substituting  and  Solving  for  D: 


|rJD  - 


4  x 


/     \ 
1 

3 


D 

=    1 


D  =  -j^-  %  1.2618 
log3 

Which  is  equal  to  the  Hausdorff-Besicovitch  dimension  of  the  Koch  curve12. 

This  is  not  a  proof  of  a  general  equation  for  the  fractal  (Hausdorff) 
dimension  of  a  self  similar  fractal  set  but  implies  that  a  general  dimension 
generating  function  is  possible: 


G(D)  =    £  (rmj 


m=l  *        ' 
Where   N  =   number    of    sides   of    the    generator. 

Where    Rm  =   ratio    of    side   m  to  the    initiator. 

Mandelbrot  claims  that  experimental  evidence  suggests  that  this  equation  holds 
whenever  this  functional  method  is  used. 

When  each  segment  of  the  generator  is  a  fixed  ratio  to  the  initiator  (as  is 
the  case  with  the  Koch  curve)  then  the  solution  to  this  equation  is  trivial: 


D         lo8(N) 


l°g(— ) 


Functional  Characterization  of  the  Koch  Method 


A  complete  and  rigorous  inductive  definition  of  most  fractal  functions 
can  stretch  the  notational  capabilities  of  the  symbolic  aspects  of  the  inductive 
and  functional   methods.   It   is  thus  generally   impractical   to  use   these   methods 


I2The     Koch     curve     was     proven     to     have     a     Hausdorff-Besicovitch     dimension     equal 

I O  £J  4 

~  1.2618  and  a  topological  dimension  equal  to  1,    Hausdorff   Dimension  und  au>seres  \l  iss 

log3 
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except  in  a  verbose  and  non-rigorous  manner.  It  can  be  insightful  however,  to 
dissect  the  beast  (once!)  and  hopefully  gain  further  intuitive  insight. 

To  simplify  the  process,  we  define  a  Koch  half-line  as  in  the  above  fractal 
but  with  a  single  line  segment  as  the  Initiator  (versus  the  equilateral  triangle). 
We  consider  the  line  interval  of  [(0,0), (1,0)]  within  the  set  R2.  This  restricts  the 
fractal  shape  that  is  drawn  to  a  partial  function  on:  [0,1]  X  [0,1]  -*  [0,1]  X  [0,1]. 

Using  the  inductive  process  to  define  the  essence  of  the  fractal 
sequencing,  we  have  the  following  definition: 

(Basis) 

Step  0 

n0  =   {(o.o,o.o), (i.o,o.o)} 

(Induction) 

Step  k 

Label  the  ordered  set  of  2-tuple  points  from"  Step  k-1  as: 

n*-i  =  {Pk-\Pk-\Pk-\,..,Pk-\  } 

where 
n    =  4(*-l)  +  i 

Determine  the  new  set  of  ordered  points  for  Step  k  as: 

n*   =  {P*1P*2P*3,...,P*m  } 

where 

m    =  4*    +1 

where 

iJfc        _      r>k-l 


P*!    =    PK~\ 

Pk5   =    Pk'\ 
Pko  =    Pk~\ 


p  k  _     p  k  -  1 
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where 


(p\,p*7,p*8)  =  Pk~\  [2R3]p*  ' 


C*  10/*  11/*  12)  -  Pk~X: 


3R4 


k-\ 


I  pk  pk  pk  ok -1 

V-*       m-3,r      m-2,r      m  - 1 


l)=    ^ 


n-1 


n-l^nj 


*-] 


and  where  the  relation 


,-iR. 


is  the  geometrical  relationship  between  the  two 
end  points  of  the  initiator  line  segment  from  step  k-1  and  the  five  points  of  the 
generator  that  compose  the  ordered  subset  for  step  k  as  above.  For  purposes  of 
brevity,  the  full  functional  definition  for  this  geometrical  relationship  is  explained 
in  detail  in  Chapter  4. 

(Extremal) 

No  point  is  an  element  of  Q  unless  it  can  be  shown  to  be  so  in  a  finite 
number  of  applications  of  clauses  1  and  2  above. 

It  is  thus  possible  (but  tedious)  to  rigorously  characterize  the  Koch 
fractal  set  within  the  well  defined  constructs  of  induction  and  the  functional 
technique. 

6.     Another  Fractal  Set,  Cantor's  Dust 

Cantor  developed  a  function13  which  used  the  same  functional  technique 
as  the  Koch  curve  but  with  a  reverse  twist.  Cantor's  function  takes  an  initiator 
(the  unit  interval  [0,1])  and  dissolves  it  into  a  discontinuous  set  which  is  as  rich 
in  points  as  the  interval  [0,1]  but  contains  no  interval  itself  [Ref.  6:pp.  22-23]. 
The  points  of  the  set  are  all  distinct  but  the  set  has  the  same  cardinality  as  the 
unit  interval.    The  best  description  for  this  set  is  that  it  resembles  a  dust. 


is 


Also  referred  to  as  Cantor's  discontinuum  or  Cantor's  triadic  set. 
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Cantor's  Dust  is  difficult  to  demonstrate  on  a  graphics  display  because  it 
quickly  dissolves  below  the  resolution  of  the  display.  Refer  to  Figure  2.2  to 
visualize  the  initiator  and  the  generator.   The  initiator   is  a  line  interval   [a  ,6  ] 

where  0  ^  a   <  b   ^1  and  the  generator  is  two  intervals  each  —  the  length  of 


the    initiator   such   that    interval    1    is    [a  ,  a  + 


(b  -a)  x  — 


and    interval   2    is 


/ 


[6- 


(b  -a)  x  — 
V  '        3 


,6].    After  the  initial  application  of  the  generator  to  the  unit 


1  2 

interval  there  will  be  two  intervals  in  the  current  construction  [0. — ]  and  [ — ,1]. 

The     second     iteration     of    the     recursive     routine     will     yield     four     intervals 


2    1 


I     t 


8 


[0.-1,  [-,-].  [-,-]  and  f-,ll. 
1     9J    l9    3J      3    9J  l9     J 

Every  initiator  that  is  created  by  the  generator  has  an  infinite  sequence 
that  is  begun  at  the  next  application  of  the  generator.  This  causes  a  series  of 
convergent   sequences   to  each   end  point  of  each   initiator.   Thus  each   initiator 
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Figure    2.2    Cantor's    Dust 
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spawns  a  convergent  sequence  toward  its  endpoints.  But  for  each  initiator  there 
are  also  two  spawned  sequences  toward  the  center.  It  is  possible  to  imagine  this 
as  an  infinitely  dividing  organism  which  leaves  behind  four  points  (eggs)  each 
time  it  divides  and  then  each  egg  itself  replicates    •  •  ■    ad  infinitum. 

This  functional  method  of  dissolxnng  a  line  is  very  powerful  as  a  tool  in 
computer  graphics  because  it  can  be  used  to  cause  many  special  effects  from 
ordinary  objects.  Mandelbrot  has  used  variations  of  this  method  to  create  images 
of  star  systems  for  example. 

The  topological  dimension  of  Cantor's  dust  is  0  and  by  Mandelbrot's 
dimension  generating  function  we  have  the  fractal  (Hausdorff)  dimension  equal 
to: 

N    /      m 
G(D)  =     £    rj      =   1 

m=l  V        ' 

Where   N  =  2 

Where   Rm  =  - 
m  3 


Substituting  and  solving 
D 


2  x 


'i 

3 

v     / 


=  1 


log  2 
D  =  -^£_=.  «  .6309 
log  3 

Cantor's    dust    was    proven    to    have    a    Hausdorff-Besicovitch    dimension    of 

-i^Jl  ~    6309  [Ref.  l:pp.  77-78]. 
log  3 

7.      A  Note  on  the  Concept  of  Similarity  and  Fractals  Sets 

The  use  of  the  mathematical  concept  of  similarity  continually  shows  up 
in  the  investigation  of  fractal  functions.  This  can  be  expected  because  of  the  type 
of  functional  building  tools  that  are  used  for  most  of  these  sets.  The  relationship 
of  the  generator  to  the  initiator  has  similarity  built  in.  Thus  Mandelbrot  is  able 
to  use  similarity  properties  in  formulating  a  dimension  generating  function  and  in 
making  claims  about  the  set's  inherent  structure. 
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The  functional  method  invented  by  Koch  and  Cantor  is  but  one  method 
of  determining  a  fractal  set.  Similarity  fractals  may  eventually  be  grouped  into  a 
class  of  fractals  (important  but  restricted).  This  type  of  functional  method 
provides  a  vehicle  for  the  creation  of  disciplined  fractal  sets  but  makes  you 
wonder   about   the   rest   of  the   space  that   fills   the  gaps  between  the  standard 

Euclidean  shapes14,  the  self-similar  fractals  and  the  infinitude  of  R1  . 


Much  like  the  transcendental  numbers  (n  and  e  for  example  )  in  all  their  uncountably  rich 
expanse  fill  the  gaps  between  the  familiar  (and  well  behaved)  sets  in  R. 
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III.    THE  IMPLEMENTATION  OF  FRACTALS  IN  COMPUTER  GRAPHICS 

This  chapter  introduces  the  reader  to  the  practical  aspects  of  implementing 
Fractals  on  a  two  dimensional  computer  graphics  display.  We  are  forced  to 
confront  the  issues  of  economy  of  scale  between  the  infinite  fractal  function  and 
the  finite  computational  environment  of  the  computer.  L'nderstanding  these 
compromises  is  a  necessary  bridge  to  successful  fractal  programming. 

A.    THE  IMPLEMENTATION  PROBLEM 

1.     Infinite  Recursion,  Stacks  and  Data  sets 

It  is  naive  to  view  the  computer  as  a  truly  infinite  abstract  machine 
which  is  capable  of  any  binary  computation  of  any  length.  Infinity  is  a  concept 
that  when  applied  to  physical  objects  quickly  breaks  down  as  soon  as  that  object 
is  bounded  in  any  way. 

It  is  possible  to  model  infinite  behavior  in  computers  though  mathematics 
(automata  theory)  and  gain  useful  insight  into  possible  capabilities  of  the 
computer  (the  use  of  push  down  automata  in  compilers  for  instance).  But  in 
order  for  automata  theory  to  make  assertions  it  is  often  necessary  to  make  the 
assumption  that  the  automata  (computer  or  an  abstract  machine)  is  in  fact 
infinite.  When  the  assumption  of  infinity  is  made,  there  are  many  powerful 
mathematical  tools  which  can  be  brought  to  bear  upon  non-intuitive  abstract 
problems  that  would  otherwise  be  functionally13  intractable  if  the  automata  was 
considered  to  have  an  arbitrarily  large  but  bounded  space.  The  question  arises, 
can  we  consider  the  implementation  of  these  constructs  (insights)  as  valid?  The 
answer  is  yes.  but  only  in  the  context  of  some  bounded  space  (for  instance,  the 
maximum  stack  space  for  a  push  down  compiler). 

The  question  of  how  many  valid  programs  can  be  recognized  by  such  a 
compiler  is  usually  too  difficult  to  determine.  The  compiler's  solution  to  the 
problem  is  to  use  a  passive  sensor  to  detect  stack  overflow  and  notify  the  user 


15  Functionally  in  the  sense  that  the  problem  may  be  decidable  but  the  solution  space  is  so 
large  and  undefined  that  it  could  not  be  determined  in  a  reasonable  amount  of  time. 
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that  his  problem  is  to  large  for  the  current  stack  space.  Although  the  compiler  is 
theoretically  a  recognizer  for  an  infinite  set  of  programs,  the  finiteness  of  the 
computer  is  the  grounding  factor. 

A  similar  paradigm  exists  for  implementing  infinite  functions  like  fractals 
in  computers.  Fractal  functions  use  recursion,  randomness,  massive  floating  point 
computations  and  large  amounts  of  primary  memory.  In  order  to  use  the  fractal 
function  productively,  we  must  manage  the  methods  we  use  and  produce  a  finite 
approximation  of  the  fractal  set  we  create.  This  can  be  done  by  a  passive  or 
active  means. 

2.     The  Bounded  Stack  Limits  Recursion 

The  Koch-like  fractal  functional  method  is  by  definition  a  recurrence 
equation.  All  fractal  methods  introduced  in  this  thesis  have  a  similar  recurrence 
definition.  Thus  it  is  impossible  to  avoid  the  use  of  recursion  in  the  production  of 
fractal  graphics16. 

Recursion  on  most  Von  Neuman  type  computers  is  implemented  on  the 
system  stack  (which  may  be  hardware  (fixed)  or  software  (in  primary  memory 
and  therefore  expandable)).  Such  a  stack  is  always  bounded  by  a  fixed  upper 
limit  of  allocatable  memory  space.  The  formal  definition  of  the  Koch-like  fractal 
method  has  no  recursive  terminating  event.  In  order  to  use  these  methods,  we 
must  determine  the  precision  that  we  require  and  develop  a  termination  event  to 
signal  the  beginning  of  the  recursive  ascent  when  this  precision  is  met. 

Failure  to  manage  the  recursive  descent  inevitably  causes  the  exhaustion 

of  the  system  stack  which  stores  the  program  state  on  each  recursive  call17.  The 
programmer  should  be  careful  to  keep  his  local  variables  at  a  minimum  in  the 
recursive  subroutine.  By  doing  so,  the  total  data  stored  during  each  recursive 
call  is  reduced. 


16  Why  should  you  try  ?  The  recursive  part  of  fractals  functions  is  the  mathematically  beau- 
tiful aspect  which  makes  them  so  eloquent  and  powerful. 

This  could  be  the  recursive  termination  event  if  you  are  not  careful.  It  might  be  a  good 
idea  to  try  this  experimentally  on  your  computer  to  determine  this  maximum  recursive  descent 
distance. 
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3.     The  Computer  Graphics  Set  Paradigm 

The  paradigm  that  is  used  for  displaying  objects  on  the  raster  graphics 
screen  is  inherently  discrete  and  two-dimensional.  A  typical  system  consists  of 
primitives  that  allow  the  user  a  view  of  his  modeling  world  as  a  largely 
unrestricted  three  dimensional  space.  These  primitives  are  limited  in  their  power 
and  are  usually  based  on  the  Euclidean  geometry  (lines,  points,  polygons  etc.). 
The  user  is  required  to  supply  viewing  parameters  that  define  a  limited  three 
dimensional  viewing  space.  These  commands  are  processed  by  the  graphics 
system  which  projects  the  objects  contained  in  the  viewing  space  onto  a  two- 
dimensional  space  that  is  the  display  screen.  The  display  screen  can  be  divided 
into  discrete  entities  called  pixels,  each  of  which  represents  one  point  on  the 
screen.  The  number  of  pixels  defines  the  resolution  of  the  screen. 

Given  the  graphics  paradigm,  it  is  natural  to  model  objects  that  are 
defined  on  the  display  as  a  collection  of  discrete  points.  When  viewed  on  the 
screen,  these  pictures  can  accurately  model  objects  of  considerable  complexity 
[Ref.  l].  This  view  is  a  departure  from  the  normal  view  of  the  graphics 
application  programmer.  The  application  programmer  views  the  objects  he 
creates  through  a  collection  of  Euclidean  geometric  primitives  (lines  polygons 
etc.)  that  are  abstracted  from  the  actual  display18.  Both  views  are  useful  models 
in  describing  objects  but  for  the  fractal  graphics  programmer  the  former  is  the 
more  powerful  because  it  neatly  maps  to  the  fractal  method  of  object 
construction. 

For  this  study,  we  view  the  world  and  the  objects  it  contains  as 
collections  (sets)  of  discrete  three  dimensional  points.  When  an  object  is  built 
through  fractal  techniques,  the  length  of  the  pixel  is  a  natural  termination  point 
for  the  fractal  recursive  process.  This  yields  an  attractive  bijective  mapping  from 
our  object  to  the  display  screen. 


18  All  computer  graphics  is  an  abstraction;  realism  is  achieved  through  deception  of  the  eye 
by  a  very  small  collection  of  colored  points.  Our  goal  then,  is  to  find  the  most  efficient  method  of 
defining  those  points. 

37 


4.     Summary 

From  our  previous  discussion,  it  should  be  obvious  that  fractal  sets  do 
not  exist  in  our  computer.  The  sets  that  we  create  are  finite  approximations  of 
the  actual  set.  The  same  is  true  of  the  fractal  pictures  that  we  create  and  to 
which  we  ascribe  fractal  dimensions.  It  is  an  illusion  of  the  eye  that  we  create  by 
taking  the  fractal  computations  below  the  resolution  of  the  screen. 

B.    MAPPING  FRACTALS  TO  A  BOUNDED  SPACE 

In  view  of  the  graphics  set  paradigm,  it  behooves  us  to  develop  a  completely 
bounded  set  space  which  has  a  one  to  one  correspondence  between  the  points 
computed  by  the  fractal  equations  and  the  entities  of  the  screen.  This  abstraction 
is  appropriate  because  these  entities  called  pixels  (which  are  nothing  more  than 
colored  points)  are  the  only  components  of  the  picture19.  The  only  reason  for  not 
using  this  methodology  before  is  that  it  was  functionally  difficult  to  compute 
(without  some  simplifying  abstractions)  the  large  number  of  pixels  in  the  display 
set. 

1.     Fractal  Recursive  Termination  Event 

Most  graphics  software  packages  provide  a  means  by  which  the  user  can 
define  the  space  on  the  display  screen  onto  which  he  wishes  to  map.  This  is  done 
by  providing  to  the  graphics  system,  viewing  parameters,  which  define  a  bounded 
three  dimensional  space  which  is  then  projected  onto  the  two  dimensional  screen. 

Many  mappings  are  possible  within  the  above  constructs  but  for  the 
fractal  programmer  it  is  convenient  to  establish  a  mapping  that  provides  for  a 
one  to  one  (or  bijective)  relationship  between  the  3  pixel  set  and  the  real-valued 
coordinate  space  in  which  we  compute  a  fractal  object.  We  develop  this  bijective 
relationship  by  defining  a  fractal  part  to  be  determined  when  the  distance 
between  generating  points  (in  real  coordinates)  is  less  than  the  distance  of  one 
pixel  in  the  mapping  of  real  to  screen  coordinates  (SCR),  Figure  3.1.  This 
mapping  is  explicit  when  the  window,  viewport  and  z- clipping  are  defined.  This 
mapping  limits  the  number  of  pixels  that  can  be  associated  with  a  given  fractal 
initiator  to  a  one  to  one  mapping  of  fractal  parts  to  the  pixel  set.    The  size  of  the 


In  fact,  this  is  the  way  the  Euclidean  graphics  abstracts  such  as  a  line  are  mapped  to  the 
screen.  At  some  point  in  the  display  process  this  mapping  must  take  place. 
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pixel    is   a  cubic   space  which   clumps  together   the    infinite   fractal   space   into   a 
discrete  number  of  cubic  spaces  equal  to  the  number  of  pixels  in  the  pixel  set. 

The  fractal  programmer  should  choose  his  viewport  and  window  carefully 
so  as  to  elicit  the  most  attractive  mapping  possible  between  his  object  and  the 
screen.  The  viewport  (SCR)  rectangle  should  be  defined  such  that  the  ratio  of 
the  X  distance  (horizontal)  to  the  Y  distance  (vertical)  is  equal  to  the  ratio  of  the 
X-Y  distances  of  his  (real-valued)  window.  The  Z-distance  (front  and  back 
clipping  plane  distance)  should  be  equal  to  the  distance  in  the  window 
coordinates  that  defines  the  length  of  the  pixel  multiplied  by  the  desired  Z-depth. 
Figure  3.2.  Formally: 
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If  the  above  relationship  holds  then  the  size  of  the  pixel  in  world  coordinates 
equal  to  the  X  (or  Y)  world  distance  divided  by  the  number  of  pixels,  Figure  3.2. 

PIX  SIZE  = 


XsCR 

The  front  and  back  Z-clipping  planes  should  be  established  such  that  the  distance 
from  front  to  back  is  equal  to  the  desired  Z  depth  (expressed  in  pixels)  multiplied 
by  the  pixel  size  (in  world  coordinates),  Figure  3.2. 

If  the  above  viewing  relationships  hold,  then  the  inclusion  of  a  simple 
check: 

IF  (InitiatorJDistance  <  PIX_S1ZE)  THEN 

return; 
ELSE 

continue; 

terminates  the  current  recursive  descent  and  begins  backtracking. 
2.     Memory  Requirements 

The   amount   of  main   memory   required   by   the   fractal   programmer   is 
directly  proportional  to  the  size  of  the  pixel  space  that  is  being  used  to  display  a 
fractal  figure  and  the  sophistication  of  the  desired  display, 
a.      Data  Locality  During  Recursive  Descent 

If  the  fractal  in  question  is  a  1-2  dimensional  display  that  is 
displayed  as  a  collection  of  points  on  a  two  dimensional  plane  then  the 
requirement  to  store  the  3-tuple  points  is  eliminated  ((x  .y  ,z  )  coordinates  which 
represent  a  point  in  R3).  The  points  can  be  computed  by  the  function  and 
displayed  on  the  screen  in  the  immediate  mode  and  then  destroyed.  The  only 
requirement  for  memory  is  the  data  that  is  germane  to  the  program  ("globals" 
and  thus  fixed  at  run  time)  and  the  amount  of  data  pushed  onto  the  run-time 
stack  for  the  recursive  calls  up  to  the  point  of  the  deepest  recursive  level.  This 
requirement  is  minimal  and  does  not  present  any  significant  limitations. 

An  example  of  such  a  fractal  is  the  Koch  curve.  The  algorithm 
presented  in  chapter  4  uses  the  data  that  is  local  to  the  subroutine  that  is 
computing  the  current  generator  with  the  input  coordinates  of  the  initiator.    The 
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inductive  nature  of  the  fractal  method  provides  data  independence  so  that  the 
programmer  does  not  have  to  have  available  to  his  subroutine  all  points 
computed  thus  far. 

b.     Memory  Requirements  for  the  Fractal  Set  Paradigm 

If  the  programmer  requires  a  more  sophisticated  display  that  includes 
a  requirement  for  hidden  surface  removal  and  (or)  lighting  enhancement  then  the 
fractal  method  becomes  just  a  part  of  the  overall  display  process.  The 
programmer  has  the  fractal  process  compute  the  points  of  the  picture  and  these 
points  are  stored  in  a  memory  structure.  This  structure  is  processed  by  a  hidden 
surface  algorithm  and  a  lighting  algorithm  and  is  then  projected  onto  the  screen. 

The  ideal  Fractal  computer  would  have  an  exhaustive  memory  and 
unlimited  computing  power  to  be  able  to  allow  us  to  store  the  entire  pixel  set  in 
memory.  This  however,  can  stretch  the  capabilities  of  most  present  day 
computers.  For  example;  if  the  pixel  set  is  defined  such  that  its  dimensions  are 
1000x1000x1000  (which  is  certainly  reasonable),  the  amount  of  storage  required 

would  be  an  array  of  109  x  24  bits20  or  24  billion  bits. 

Current  techniques  exist  where  this  storage  can  be  minimized.  A  Z- 
buffer  array  can  be  used  to  store  the  current  Z  coordinate  of  the  forward  most 
displayed  pixel.  By  using  this  method  the  fractal  computation  can  be  made  and 
then  a  hidden  surface  computation  can  be  immediately  invoked  to  check  the 
visibility  of  the  point  by  checking  it  against  any  other  point  in  the  same  position 
on    the    X-Y    plane.    This    technique    reduces    the    above    space    requirements    to 

106  x  24  bits  +  (Z- buffer  =  106  x  32  bits)21  or  56  million  bits. 

A  fractal  programmer  must  manage  his  memory  resources  carefully 
in  order  to  maximize  the  computational  resources  available  to  him.  This  usually 
requires  a  tradeoff  between  efficiency,  realism  and  the  memory  space  utilized. 
The  algorithms  for  hidden  surface  and  lighting  are  well  covered  in  the  literature 
and  although  they  are  integral  aspects  of  the  fractal  realism  issue  they  are  not 
germane   to    this    study.    Throughout    the    remainder   of  this   study,   we   are   not 


"    24  bits  for  a  machine  assuming  the  RGB  color  system  with  24  bit  planes.    Each  color;  red, 
green  and  blue  would  then  have  8  bits  of  precision. 

1  Assuming  the  Z  coordinate  is  stored  as  a  32  bit  floating  point  number:  this  number  can  be 
further  reduced  if  we  restrict  the  Z  precision. 
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main() 
LOAD  INITIATOR  AND  GENERATOR  DATA 


concurrent; 
generate  (INITIATOR  l) 
generate  (INITIATOR  2) 
generate  (INITIATOR  3) 


generate(INITIATOR  n); 
end  concurrent; 
end  main; 

generate(INITIATOR  I) 

BUILD  GENERATOR  FROM  INITIATOR 

if  (Distance  <  PIX_size)  then 
PLOT  POINT 

return; 
endif 


concurrent; 
generate  (INITIATOR  1) 
generate(INITIATOR  2) 
generate(INITIATOR  3) 


generate(INITIATOR  m); 
end  concurrent; 
end  generate; 


Figure  3.3.  An  Algorithm  Which  Utilizes  Concurrent  Operations. 
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IV.    FRACTAL  COMPUTATION  IN  R2 

The  algorithm  introduced  in  the  following  section  is  capable  of  computing  a 
very  broad  class  of  Koch-like  fractal  curves  within  R2.  It  provides  the  fractal 
graphics  programmer  with  the  basic  tools  for  fractal  computation  and  an 
algorithmic  template  that  can  be  used  for  many  applications.  The  graphics 
programmer  needs  to  fully  understand  fractal  programming  within  R"  before  he 
attempts  the  more  complicated  concepts  of  fractal  terrain  modeling  in  R3. 

A.    THE  GEOMETRY  OF  INITIATOR  -»  GENERATOR 

The  geometric  relationship  between  the  Koch  curve  initiator  and  generator 
can  be  described  through  a  very  simple  set  of  data  that  captures  the  essence  of 
the  Koch  curve.  The  method  introduced  also  allows  us  to  vary  the  data  which 

defines  the  generator  and  compute  many  different  fractal  shapes    . 

The  general  strategy  for  computing  generator  points  from  a  set  of  initiator 
points  is  to  determine  two  lines  that  intersect  at  an  unknown  generator  point. 
Figure  4.1. a.  The  unknown  generator  point  is  defined  by  constant  relationships 
between  the  initiator  and  generator.  The  first  line  is  the  perpendicular  from  the 
generator  point  to  the  initiator  line.  The  second  line  is  formed  between  the  first 
point  of  the  initiator  and  the  unknown  generator  point.  All  data  to  compute  the 
line  equations  are  derivable  from  the  two  endpoints  of  the  initiator  or  from 
constant  initiator/generator  ratios.  For  simplicity's  sake,  we  ignore  the  divide- by- 
zero  problems  that  are  encountered  when  any  of  the  lines  are  parallel  to  the  X  or 
Y  axis.  These  situations  only  simplify  the  computations  and  their  solution  is 
demonstrated  in  Appendix  A. 

If  we  label  the  endpoints  of  the  initiator  as  P1=(X1,Yj)  and  P:=(X:.Y:) 
(Figure  4.1. a)  then  the  slope  of  the  initiator  line  is: 

Y,  -  Y, 


Slope,  init  = 


X2-  X, 


We  must  be  careful  in  our  terminology  because  the  relationship  described  can  draw  shapes 
that  do  not  avoid  self  intersection  and  thus  must  be  considered  quasi-fract al. 
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The  slope  of  the  perpendicular  intercept   line  is  the  negative   inverse  of  the 
slope  of  the  initiator,  hence  the  slope  of  this  line  is: 

Slope,  perp  = 


Slope. init 


The  intercept  point  on  the  initiator  can  be  determined  by  using  a  constant 
distance  ratio  as  shown  by  the  following  equations: 

Xj  +  (Generator. ratio. constant    x  X2) 
X.perp  =  — 


Y.perp  = 


1  +  Generator. ratio,  constant 

Fj  +  [Generator. ratio. constant    x   Y 2) 
1  +  Generator. ratio. constant 


The  value  of  "  Generator. ratio. constant"  is  a  fixed  constant  (although  it  might  be 
interesting  to  randomize  it)  where  you  determine  at  what  point  the  generator 
intercept  point  intersects  the  initiator  and  express  it  as  demonstrated  in  Figure 
4.1.b.  This  ratio  can  be  determined  graphically,  through  hand  calculation  or  via 
an  interactive  automated  means. 

With  the  slope  and  a  point  of  the  line  (X.perp, Y.perp)  we  can  determine  the 
line  equation  for  the  perpendicular  line  by  determining  the  Y  intercept: 

Y. intercept. perp  =  Y.perp  -  (Slope. perp  x  X.perp) 

This  yields  the  line  equation  for  the  perpendicular  line: 

Y  =   (Slope. perp  x  X)  +  Y. intercept. perp 

To  determine  the  generator  line  equation  for  the  line  segment  which  connects 
the  first  point  of  the  initiator  and  the  unknown  generator  point,  we  need 
constant  information  about  the  angle  between  the  initiator  and  this  line.  This 
angle  is  always  constant  with  respect  to  the  initiator  and  like  the  ratio 
information  above,  it  is  recorded  as  a  constant  at  run  time  (see  Figure  4.1.c), 
note  that  the  angle  may  be  positive  or  negative. 
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Fig  4.1. a  Intersecting  Lines  Determine  a  Generator  Point. 
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Fig  4.1.b.  The  Generator  Ratio  Constant 
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Fig  4.  I.e.  The  Angle  0   '\/*a 


Figure  4.1  Building  the  Generator  with  Intersecting  Lines 
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We  record  the  data  about  the  angle  0  as  the  tangent  of  the  angle.  With  this 
information,  the  slope  of  the  generator  line  can  be  determined  with  the  following 
equation: 

01  Tan  6    +  Slope. init 

Slope. gen  = 


(1  —   Tan  6  )  x  Slope. init 
The  Y  intercept  for  the  generator  line  can  now  be  determined: 

Y. intercept. gen  =  Yl  —  (Slope. gen  x  Xx) 

This  yields  the  line  equation  for  the  generator  line: 

Y  =  (Slope. gen  x  X)  +  Y. intercept. gen 

The  Cartesian  points  of  the  unknown  generator  point  can  now  be  determined 
by  intersecting  the  two  line  equations  and  solving  for  X  (Figure  4.1. a)  then 
substituting    X  into    one    of  the    line    equations    and    solving   for   Y       .    The 

equations  follow: 

v  ,.  .  Y. intercept. gen 

A        =   Y. intercept. perp  —  

8  Slope. gen  —  Slope. perp 

Ygen  =  (Slope. gen  x  Xgen)  +  Y. intercept. gen 

The  constant  data  for  the  Koch  curve  that  corresponds  to  this  geometric 
method  is  illustrated  in  Figure  4.2. 

There  are  many  different  ways  to  build  a  generator  given  an  initiator  using 
standard  geometric  constructs  but  the  method  introduced  allows  experimentation 
with  the  Koch  function  to  discover  new  shapes.  By  varying  the  tangent  of  6  and 
the  fixed  ratio  of  Figure  4.1  we  can  describe  new  generator  constructions.  These 
new  constructions  can  be  used  in  the  same  algorithms  that  compute  the  Koch 
curve.  We  initiate  the  recursion  on  a  generating  structure  built  of  line  segment 
initiators  and  allow  the  recursion  to  progress  until  a  terminating  event25,  creating 
many  diverse  shapes.  Figures  4.6  through  4.9  demonstrate  some  of  these  shapes. 


The  length  of  a  pixel  or  an  arbitrary  line  length. 
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B.    THE  MID-POINT  DISPLACEMENT  TECHNIQUE 

Mandelbrot  [Ref.  l:pp.  43  and  pp.  233-234]  uses  an  alternate  technique  to 
draw  the  Koch-like  curves  that  is  equally  valid.  He  calls  his  technique  mid-point 
displacement  because  he  determines  the  generator  via  fixed  relationships  that 
displace  a  point  from  the  mid-point  of  the  initiator.  This  method  uses  many  of 
the  same  geometric  relationships  that  are  used  above  but  provides  a  different 
progression  to  the  method  of  building  the  generator.  This  new  view  allows  us  to 
look  at  the  relationship  between  the  initiator  and  generator  in  a  slightly  different 
light.  By  so  doing  we  are  provided  new  insight  as  to  how  we  might  alter  the 
relationship  to  create  new  images  and  sets. 

The  mid-point  displacement  technique  can  be  useful  for  two  other  reasons.  It 
is  the  best  known  method  for  fractal  set  building  and  because  of  this,  it  facilitates 
communication  between  fractal  programmers.  Most  of  the  terrain  models  that 
have  been  developed  use  the  mid-point  technique.  It  also  provides  an  easier 
method  to  avoid  line  intersection. 
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The    Midpoint   Displacement 
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Figure  4. 3. a  First  Midpoint  Displacement 


Figure  4.3.b  Second  Midpoint  Displacement 


THE  KOCH  GENERATOR 


V--- 


^\ 


Figure  4.3.c  Third  Midpoint  Displacement. 


Figure  4.3.  The  Koch  Generator  Using  Midpoint  Displacement. 

The  mid-point  displacement  method  is  demonstrated  in  Figure  4.3.  The 
method  progresses  by  taking  the  initial  initiator  and  applying  the  first  midpoint 
displacement.  This  yields  the  figure  demonstrated  in  Figure  4. 3. a.  The  next  step 
performs  a  mid- point  displacement  on  the  left  initiator  created  by  step  1.  This 
yields  the  figure  demonstrated  in  Figure  4.3.b.  The  third  and  final  step  is  to 
replace  the  right  initiator  created  by  step  1  as  demonstrated  in  Figure  4.3.c. 
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The  inversion  of  the  direction  of  the  mid-point  application  (in  Figure  4. 3. a 
the  displacement  is  above  the  initial  initiator  where  in  Figure  4.3.b  it  is  below) 
can  be  accomplished  with  a  single  computational  procedure.  We  only  need  to 
invert  the  position  of  point  1  and  point  2  in  the  parameters  of  that  procedure. 
The  parameter  inversion  changes  the  orientation  of  the  initiator  in  space  with 
respect  to  the  computation  of  the  midpoint  displacement.  The  procedure  blindly 
computes  the  mid-point  displacement  relative  to  a  fixed  relationship  to  the 
initiator  input  points.  This  operation  implicitly  defines  the  orientation  of  the 
generator  in  space. 

The  geometry  for  computing  the  midpoint  given  two  initiator  points  can  be 
computed  in  precisely  the  same  manner  as  the  intersecting  line  algorithm  (The 
Koch  midpoint  ratios  are  1.0  and  the  angle  0  is  .437  radians).  An  alternative 
method  is  to  use  the  equations  of  a  line  (or  a  plane)  normal.  The  second  method 
(utilizing  the  normal)  provides  a  geometric  relationship  which  is  intuitively 
appealing.  Its  appeal  comes  from  the  desire  to  modify  the  length  of  the 
displacement  relative  to  the  initiator  with  a  random  scaling  factor. 

The  mid-point  displacement  technique  has  some  advantages  over  the  line 
intersection  algorithm.  Random  modification  of  the  length  of  the  displacement 
along  the  normal  (from  the  computed  generator  point  to  the  initiator  using  the 
line  intersection  algorithm)  is  not  intuitively  appealing.  It  requires  a  translation 
of  the  desired  displacement  into  an  angle  (the  angle  9  (Figure  4.1.c)  between  the 
initiator  line  and  the  unknown  generator  intercept  line).  The  control  of  that 
angle  is  less  intuitive  than  the  control  of  a  displacement  length.  The  geometry 
for  mid-point  displacement  using  the  initiator  normal  is  introduced  in  Chapter  5. 

C.    A  KOCH-LIKE  FRACTAL  ALGORITHM 

Implementing    the    above    function    is    a    relatively    easy    process    that     is 

demonstrated  in  this  section  and  Appendix  A  via  a  gradual  unwrapping  of  the 
layers  of  complexity  that  are  required  to  successfully  implement  the  algorithm. 
The  algorithm  roughly  follows  the  template  used  in  Chapter  3  to  demonstrate 
concurrent  processing.  A  C-like  language  is  used  for  the  algorithms. 

The  first  algorithm  (Figure  4.4)  is  a  template  that  delineates  the  basic 
processing  steps.  This  recursive  process  is  typical  of  fractal  functions  and  can  be 
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used  as  a  template  for  many  fractal  programs.  The  second  algorithm  (Figure  4.5) 
is  an  expansion  of  the  first  and  demonstrates  the  replacement  of  a  given  initiator 
using  the  line  intersection  method. 

Appendix  A  is  a  complete  Fractal  program.  This  program  was  used  to 
produce  the  data  for  Figure  2.1  and  Figures  4.6  through  4.9.  This  program 
demonstrates  the  precautions  that  must  be  taken  to  avoid  divide- by- zero  when 
lines  are  parallel  to  the  X  or  Y  axis. 


main() 

{ 

Load  Initiator  Coordinates; 

Load  Generator  Relationship  Values; 

for  1=1;  I<=  Number  of  Initiators;  1=1+1; 

{ 
generate(X(I)1,Y(I)1,X(I)2,Y(I)2) 

} 
} 

generate(X1,Y1,X2,Y2) 
{ 

Determine  Distance  Between  the  Endpoints  of  the  Initiator 

if   (DIST  <  Pixel.length)    return; 

Replace  Initiator  with  the  Computed  Generator; 

Load  the  New  Initiator  Data  into  Local  Generator  Array 

for  J=l;  J<=  Number  of  Generator  Segments;  J=J+1; 

{ 
generate(Xgen(J)!.  Ygenp^,  Xgen(J)2,  Ygen(J)2); 

} 
} 

Figure  4.4.  High-level  View  of  the  Koch  Algorithm. 
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main() 

( 

Load  Initiator  Coordinates; 

Load  Generator  Relationship  Values; 

for  1=1;  I<=  Number  of  Initiators;  1=1+1; 
{ 

} 
} 

generate(X1,Y1,X2,Y2) 
{ 

/*  Determine    Distance   Between   the    Endpoints   of    the    Initiator    '  / 
DIST  =  sqrt((X2-X1)**2  +  (Y2-Y1)**2); 

/*  //    Distance   is   less    than   Pixel   Length;  Plot    and  Return    */ 
if   (DIST  <  Pixel. length) 

{ 
plot.point();        /*    Your    Graphics    Point   Plotting   Routine    */ 

return;  /*  Point   1  and  2  plot  the   same   pixel   */ 

> 

I*  Load    The   Endpoints   of    the    Initiator    into   the    Generator    Array    '/ 
Generator. X[0]  =  Xp 
Generator.  Y[0]  =  Y,; 

Generator. X[Number.of.generator. points  +  l]  =  X:; 
Generator. Y[Number.of.generator. points  +  l]  =  Y2; 

/*  Determine    Slope    of    the    Initiator    */ 
Slope. init  =  (Y2-Y,)  /    (Xj-XJ; 

-  continued- 


Figure  4.5.  Detailed  \'iew  of  the  Koch  Algorithm. 
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/*  Calculate    the    Unknown   Generator   Point   via   Intersecting  Lines   */ 
for  J=l;  J<=  Number.of.generator. points;  J=J+1; 

{ 

/*  Determine   the    Generator   Point   Intercept   on   the   Initiator   */ 
X.perp  =   (Xj  +  (G enerator .ratio. constant  [J]  *  X2))  / 

(1+  Generator. ratio. constant  [J]); 
Y.perp  =  (Yl  + [Generator. ratio. constant  [J]  *  Y2))  / 

(1+  Generator. ratio. constant  [J]); 

/*  Determine    the    Slope    of    the   Perpendicular    Line    */ 
Slope. perp  =  (-1  /  Slope. init); 

/*  Determine    the    Y -intercept    of    the    Perpendicular    Line    */ 
Y. intercept. perp  =  Y.perp  -  (Slope. perp  *  X.perp); 

/*  Determine    the    Slope    of    the    Generator   Line    */ 
Slope. gen  =  Generator. tan. theda  [J]  +  Slope. init)  / 
(1- Generator. tan. theda  [J]  *  Slope. init); 

/*  Determine    the    Y  -  intercept   of    the    Generator   Line    */ 
Y. intercept. gen  =  Yl  -  (Slope. gen  *  X2); 

/*  Determine    the    Unknown    Generator   Point   */ 
Generator. X[J]  =  (Y. intercept. perp  -  Y. intercept. gen)  / 

(Slope. gen  -  Slope. perp); 
Generator. Y[J]  =  Slope. gen  *  Generator. X[J]  +  Y. intercept. gen; 

} 

for  K=0;  K<=  Number.of.generator. points  +  1;  K  =  K  +  1; 

{ 
generate(G  enerator. X[K], Generator.  Y[K], 

Generator. X[K+l],Generator.Y[K  +  l]); 

} 
}     /*  End   Generate    */ 


Figure  4.5.  Detailed  View  of  the  Koch  Algorithm  (continued). 
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D.  IMPLEMENTATION  STRATEGIES 

There  are  numerous  ways  to  display  the  fractal  shapes  that  the  above 
algorithm  is  capable  of  computing.  The  graphics  primitives  required  are  limited 
to  the  standard  initiation  and  termination  commands  coupled  with  the  ability  to 
plot  a  point  (or  alternatively  a  line).  Any  raster  graphics  system,  plotter  or 
similar  technology  suffices. 

The  algorithm  can  be  extended  to  include: 

-  Online  generator  drawing  to  compose  a  generator  relationship  visually. 

-  Rotation  in  3  dimensions  (if  your  system  has  this  capability). 

-  Variation  of  the  inductive  application  of  the  generator  by  the  inclusion  of 
randomness  with  respect  to  the  generator  constants. 

Figures  4.6  through  4.9  represent  a  few  of  the  shapes  that  this  algorithm  can 
compute.  Each  figure  has  the  generator  data  used  to  compute  the  shape  and  a 
progression  that  shows  the  first  two  recursive  iterations. 

E.  SUMMARY 

The  line  intersection  algorithm  as  it  stands  is  not  very  useful  for  the 
production  of  graphics  images  of  realistically  textured  terrain.  Its  importance 
results  from  its  encapsulation  of  the  essence  of  the  non-random  inductive  fractal 
method.  This  algorithm  demonstrates  the  idea  and  intent  of  fractal  functions  and 
their  implementation  within  computer  graphics.  The  potential  fractal 
programmer  must  throughly  understand  the  salient  parts  of  this  chapter  before 
successfully  attempting  fractal  images  in  three  dimensions  (i.e.  before  climbing  the 
mountain  Chapter  5). 
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A   Cloud-like    Shape 
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1st    Iter&ti 


2nd   Itera.ti<. 


Figure  4.6.  A  Cloud-like  Shape 
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Boxes    Ad    Infinitum 


Initiating    Structure : 
Equilateral     Triangle 


1st    Tterati 


GENERATOR  DATA 

Number    of   points    =    3 

Pointl : 
Angle    =    O.OOOO    rads 
Ratio    =    1 . OOOO 

Point2: 

Angl e    =    O. 7854    rads 
Ratio    =    l.OOOO 

Point3: 
Angle    =    0.4636    rads 
Ratio    =    200.00 


Figure  4.7.  Boxes  Ad  Infinitum 
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An   Exaggerated   Koch   Curve 
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Figure  4.8.  An  Exaggerated  Koch  Curve 
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A  Plane    Filling    Curve 
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-continued- 


This  plane  filling  curve  requires 
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ection of  the  generator  app 1 i  cm- 
tion.  This  reversal  is  explained 
on     the    next    page . 
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Figure  4.9.  A  Plane  Filling  Curve. 
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Figure  4.9.  A  Plane  Filling  Curve  (continued) 
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V.    FRACTAL  GEOMETRY  FOR  GRAPHICS  TERRAIN 

One  of  the  most  widely  recognized  fractal  images  found  in  the  literature  is  of 
the  mountain  scene.  This  type  of  terrain  modeling  is  perfectly  attuned  to  the 
fractal  technique.  The  reason  for  this  is  that  mountains  are  highly  irregular 
shapes,  with  a  rough  but  consistent  texture  when  viewed  from  a  distant  vantage 
point.  It  is  appropriate  then,  to  introduce  graphics  terrain  simulation  techniques 
through  this  model. 

This  chapter  describes  the  theory  and  techniques  of  simulating  mountainous 
terrain  with  computer  graphics.  It  provides  the  blueprint  for  fractal  graphics 
programming  within  R3  by  providing  general  tools  and  a  methodology  that  is 
easily  adapted  to  many  other  modeling  needs. 

A.    MODELING  MOUNTAINOUS  TERRAIN 

For  the  programmer  who  fully  understands  the  essence  of  the  method  of 
fractal  programming  introduced  in  Chapter  4,  the  movement  into  programming 
in  R3  is  not  difficult.  The  primary  differences  lie  in  the  quantum  jump  in 
computing  resources  that  are  required  and  the  requirement  to  perform  the 
generator  geometry  in  R3  versus  R2.  The  theory  and  technique  of  fractals  does 
not  change  substantially. 

Chapter  3  provided  a  rough  framework  to  begin  the  coalescence  of  fractal 
programming  into  a  workable  technique.  We  need  to  develop  a  number  of  tools 
from  that  chapter  and  use  standard  compiiter  graphics  techniques  to  manage 
those  tools.  To  this  framework,  we  add  new  fractal  functions  which  provide  the 
texture  of  realism  for  our  simulated  mountain. 

1.     The  Artist's  Model 

One  way  for  an  artist  to  build  a  physical  relief  model  of  a  mountain  is  to 
use  a  framework  to  provide  structure  to  the  model  and  a  texturizing  clay  to 
provide  realism.  The  artist  might  use  chicken  wire  on  top  of  *mall  boxes  as  the 
frame  with  modeling  clay  as  the  texturizing  element.  His  choice  of  clay  is 
predicated    by    the   type   of  look    that    he   wants    to   achieve.     The   chicken   wire 
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provides  an  inexpensive  and  disguised  method  to  quickly  build  the  mountainous 
shape  and  structure.  This  method  minimizes  the  cost  and  time  to  build  up  the 
clay. 

The  artist  continues  the  modeling  process  after  the  development  of  the 
basic  mountain  shape  to  achieve  hues  and  contrast  in  the  coloration.  He  might 
achieve  this  by  the  use  of  natural  lighting  to  cast  shadows  or  by  a  careful 
painting  of  prominent  features. 

2.     The  Fractal  Programmer's  Model 

There  is  very  little  in  science  that  is  truly  new  or  innovative.  We  borrow 
the  essence  of  the  above  idea  to  guide  us  in  developing  a  model  for  the  discrete 
computation  of  our  two  dimensional  picture  of  the  mountain.  This  section 
describes  the  process  intuitively  and  leaves  the  implementation  details  to  later 
sections. 

a.     The  Lattice  Control  Structure 

The  pixel  space  that  we  developed  in  previous  chapters  can  be 
divided  into  discrete  cubic  units  by  use  of  a  concept  from  mathematics  called  a 
lattice  (in  our  case  we  can  view  it  as  three  dimensional  graph  paper).  This  lattice 
serves  as  our  controlling  structure,  the  equivalent  of  the  chicken  wire  structure 
above.  It  is  beneficial  to  build  the  lattice  as  a  structure  with  well-formed 
relationships,  where  the  number  of  lines  evenly  divides  the  boundaries  of  the  pixel 
space  and  each  line  is  a  constant  distance  from  its  neighbors.  By  this  method,  we 
do  not  have  to  store  the  lattice  but  can  express  it  as  a  mathematical  function  of 
the  pixel  space. 

The  lattice  can  be  very  useful  in  developing  a  rough  approximation 
of  the  mountain  that  we  wish  to  model.  This  can  be  done  in  many  different  ways 
but  should  result  in  a  stick  frame  model  of  the  mountain  (a  connected  polygon 
mesh  like  that  of  Figure  5.1). 

The  frame  can  be  developed  through  an  online  graphics  interface 
that  allows  the  programmer  to  select  a  ground  level  plane  of  the  lattice  and 
provide  a  means  to  visually  select  points  for  the  rough  outline  (essentially  draw 
the  framework).    This  approach  is  useful  when  a  particular  shape  is  desired. 
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A  frame  can  also  be  developed  using  fractal  functions  to  pervert  the 
lattice  into  a  controlled  random  shape  from  a  given  plane  of  the  lattice.  This  is  a 
powerful  method  that  can  be  controlled  via  bounds  on  the  random  tools, 
heuristics  or  discrete  functional  bounding  of  the  fractal  function.  This  approach  is 
most  useful  when  a  class  of  mountain  shapes  are  required  but  no  particular 
mountain  needs  to  be  modeled,  i.e.  when  random  landscapes  suffice. 
Alternatively,  the  stick  frame  of  the  mountain  can  be  determined  via  manual 
(hand  computation)  means.  This  approach  is  tedious  and  limited  in  its  flexibility, 
and  is  not  recommended. 

b.  Surface  Texture  via  Fractal  Functions 

The  next  step  in  the  creation  of  a  mountain  is  to  provide  the 
graphics  clay  to  cover  our  stick  frame  model.  This  clay  is  a  fractal  function  which 
closes  the  polygons  of  the  frame  model  with  an  inductive  process  that  provides  a 
continuous  pixel  surface  for  the  entire  structure  of  the  mountain. 

The  initiator/generator  paradigm  is  used.  The  initial  set  of  initiators 
is  the  frame  described  above  and  the  generator  is  a  similar  geometrical  shape  that 
reduces  in  size  continually  until  it  becomes  the  size  of  a  pixel  and  is  mapped. 

After  the  stick  frame  of  the  mountain  is  developed,  this  texturizing  of 
the  surface  becomes  an  automatic  process  that  terminates  when  each  geometrical 
shape  that  makes  up  the  framework  is  reduced  to  a  continuous  set  of  pixels  in  the 
pixel  set.  At  this  point,  the  mountain  exists  in  the  pixel  set  (memory)  but  must 
be  provided  color  and  light  to  bring  it  to  life. 

c.  Hidden  Surface  Elimination 

The  pixel  set  has  the  entire  stnicture  of  the  mountain  in  memory, 
but  we  can  project  only  a  two  dimensional  image  of  one  plane  onto  the  screen. 
The  fractal  function  which  texturizes  the  surface  does  not  concern  itself  with 
local  computations  so  many  overlapping  pixels  are  mapped  to  the  pixel  set. 
There  are  two  reasons  then,  why  we  need  hidden  surface  removal  (in  this  case 
better  referred  to  as  hidden  pixel  elimination).  First  we  have  to  eliminate  the 
back  or  hidden  sides  of  the  mountain  by  projecting  only  those  pixels  which  are 
visible  along  the  axis  of  sight  to  the  perpendicular  planar  surface  of  the  display 
screen.    The   second   kind   of  hidden   pixel   removal   is   caused   by   mapped    pixels 
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which  were  covered  up  by  other  recursive  fractal  descents  either  before  or  after 
the  pixel  was  mapped. 

The  removal  of  hidden  pixels  is  greatly  facilitated  by  the  use  of  the 
concept  of  the  pixel  set.  In  standard  computer  graphics  hidden  surface 
elimination,  the  programmer  is  confronted  with  graphics  primitives  which  are 
functionally  continuous  Euclidean  shapes.  To  effectively  remove  the  hidden  parts 
of  these  shapes  is  in  general  very  tedious  and  mathematically  complicated.  Since 
the  graphics  programmer  is  shielded  from  the  primitive  -*  pixel  mapping,  he  is 

functionally  denied  access  past  the  simplifying  abstraction26  of  graphics  Euclidean 
primitives.  The  fractal  programmer  must  have  access  to  this  level  of  the  graphics 
mapping  and  thus  can  use  simple  techniques  to  determine  if  a  pixel  is  hidden  or 
visible. 

The  simplest  and  most  economical  means  available  to  provide  hidden 
pixel  elimination  is  through  the  use  of  Z-buffer  algorithms.  With  this  method,  the 
determination  of  whether  a  pixel  is  hidden  can  be  appended  to  the  pixel  set 
mapping  process.  The  Z  coordinate  of  a  pixel  that  is  to  be  mapped  is  checked 
against  the  Z  coordinate  of  the  pixel  currently  in  the  pixel  set  at  the  same  row 
and  column  of  the  three  dimensional  array  used  to  store  the  pixel  set.  If  the  pixel 
is  closer  to  the  planar  surface  of  the  display  screen,  then  the  Z  coordinate  is 
changed  to  reflect  the  position  of  the  newly  mapped  pixel. 

The  Z-buffer  approach,  while  powerful,  does  limit  the  fractal 
programmer's  flexibility.  The  axis  of  sight  toward  the  mountain  must  be 
determined  prior  to  the  fractal  recursive  process  so  that  the  determination  of  the 
line  through  the  (now  two  dimensional)  pixel  set  is  known.  The  Z-buffer  becomes 
an  adjacency  matrix  to  the  pixel  set  and  can  retain  information  about  forwardly 
displayed  pixels  only.  All  information  is  lost  about  other  pixels  that  were 
computed  in  the  fractal  process.  If  another  view  of  the  mountain  is  required  then 
the  entire  pixel  set  has  to  be  recomputed  with  a  new  axis  of  sight.  If  the  fractal 
function  uses  (non-tabular)  random  techniques  then  the  mountain  varies  with 
each  view. 


*  The  abstraction  provided  by  Euclidean  primitives  is  a  powerful  one  when  the  alternative  of 
pixel  mapping  is  considered.  Without  some  powerful  mapping  tool  (such  as  fractal  functions),  the 
pixel  level  modeling  process  is  in  general  very  difficult. 
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Most  fractal  pictures  consume  such  vast  computing  resources  that 
only  one  view  is  computed  for  a  given  picture.  As  more  requirements  for  graphics 
terrain  are  determined,  a  more  powerful  method  has  to  be  used  to  retain  all  of 
the  computed  pixels  in  the  three  dimensional  pixel  set.  This  method  requires 
that  all  pixels  be  stored  in  the  three  dimensional  array  previously  described.  The 
hidden  surface  calculation  can  then  be  performed  during  the  pixel  mapping 
operation  or  as  a  separate  calculation  that  is  performed  after  all  fractal  recursion 
has  terminated. 

As  specified  in  Chapter  3,  the  full  array  approach  requires  large 
amounts  of  memory.  This  method,  however,  allows  the  computed  fractal 
mountain  to  become  an  entity  that  can  be  manipulated  versus  an  instance  of  the 
fractal  mountain  as  above. 

Both  methods  are  viable  but  the  latter  approach  provides  more 
flexibility  for  the  programmer  whereas  the  first  approach  is  a  response  to  the 
economies  of  scale  of  data  processing.  As  new  architecttires  are  developed"'  with 
capacities  geared  toward  fractal  image  computation,  the  first  method  can  be 
eliminated. 

d.     Illuminating  the  Mountain 

If  you  stood  on  the  dark  side  of  the  moon  without  illumination,  the 
mountains  and  craters  of  the  moon  would  not  be  visible.  They  still  exist  however, 
just  as  our  imaginary  mountain  exists  in  memory.  In  order  to  visualize  them,  we 
must  illuminate  them. 

Illumination  in  computer  graphics  is  achieved  by  varying  the  light 
intensities  of  pixels  displayed  on  the  screen.  The  color  mixture  of  these  discrete 
points  determines  the  lighting  effect  that  a  viewer  perceives.  This  perception  is 
not  reality  but  another  deception  caused  by  scale  and  composition.  A  lighting 
model  then,  is  one  which  is  able  to  abstract  the  essence  of  color  from  a  real  world 


27  An  Ideal  architecture  is  one  with  a  large  main  memory  and  parallel  processing  capabilities 
with  K  processing  elements  (where  K  is  greater  than  the  maximum  recursive  descent  distance). 
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object  and  transform  that  essence  into  a  set  of  color  values  (intensities)  that 
accurately  deceive  the  human  eye  via  the  graphics  medium. 

The  literature  on  computer  graphics  contains  many  lighting  models 
with  diverse  approaches  to  the  same  problem.  Many  of  these  models  (like  those  of 
hidden   surface)    concern  themselves  with   illumination   of  continuous   Euclidean 

surfaces  and  as  such,  are  not  directly  germane  to  our  study28. 

An  object  in  space  is  a  composition  of  basic  elements.  These  elements 
interact  with  the  physics  of  light  reflection  to  create  the  spectrum  of  light  that 
our  eyes  decode.  In  a  graphics  image,  this  process  has  to  be  simulated  with 
discrete  lighting  intensity  values  for  each  pixel.  Thus,  the  illumination  of  the 
mountain  is  a  two  step  process;  the  fractal  entities  that  are  mapped  to  the  pixel 
set  must  be  provided  with  a  basic  color,  and  these  colors  have  to  be  highlighted 
and  dimmed  by  the  lighting  algorithm. 

The  basic  color  can  be  determined  during  the  pixel  mapping  event  of 
the  fractal  recursion  process  or  as  a  separate  process  prior  to  or  in  conjunction 
with  the  lighting  algorithm.  This  color  can  add  realism  to  the  picture  through 
heuristics  which  the  programmer  defines.  Most  mountains  are  composed  of 
different  types  of  rocks  and  flora  and  these  elements  change  at  different  altitudes. 
This  type  of  heuristic  combined  with  some  random  control  structure  (i.e.  to  vary 
the  snow  peak)  can  provide  for  improved  realism  (versus  making  the  whole 
mountain  brown).  The  process  of  determining  the  basic  color  must  be 
accomplished  prior  to  applying  the  lighting  algorithm  since  the  lighting  algorithm 

can  only  vary  the  intensities  of  an  existing  color29.  Developing  the  process  of 
basic  color  determination  is  best  accomplished  through  trial  and  error.  It  is  the 
artistic  aspect  of  developing  fractal  mountains. 

The  general  process  of  computer  graphics  illumination  conr-rns  itself 
with  casting  shadows  from  one  object  to  another  given  a  direction  from  an 
imaginary  light  source  and  with  highlighting  surfaces  which  are  directly  exposed 
to  the  source.    A  surface  is  highlighted  relative  to  the  angle  at  which  the  light 


"8  The  Gouraud  model  (intensity  interpolation  shading)  for  instance 

Sine* 
covered  here 


Since    so   many   diverse    color    models    exist,    the   details   of  color   representation    are   not 
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source's  rays  strike  the  surface.  This  poses  special  problems  for  fractal  surfaces 
due  to  their  discontinuity  at  every  point. 

The  process  of  illuminating  a  fractal  surface  is  best  aided  by 
divorcing  the  lighting  process  from  the  fractal  computation  process  (except  as 
noted  above).  It  is  beneficial  to  view  the  pixel  set  as  a  collection  of  pebbles  which 
have  size  and  position.  This  abstraction  allows  us  to  view  the  pixel  as  a 
continuous  space  that  can  block  light  (cast  shadows)  and  for  which  an  angle  of 
illumination  can  be  determined  (usually  in  conjunction  with  neighboring  pixels). 

A  well  formed  fractal  mountain  surface  is  completely  connected  (no 
space  between  adjacent  pixels  in  the  pixel  set).  Thus  the  surface  can  also  be 
viewed  as  a  continuous  (while  very  rough)  surface  where  reflected  light  can  be 
cast  from  or  to  adjacent  pixels. 

One  lighting  model  which  fits  the  fractal  process  is  the  Torrance- 
Sparrow  model  [Ref.  8:pp.  578-579]. 

This  model  views  an  object  as  a  collection  of  facets  which  is  each  a  perfect  reflec- 
tor (i.e.  does  not  absorb  light).  The  orientation  of  each  facet  is  given  by  the 
Gaussian  probability  distribution  function  (i.e.  the  smooth  surface  of  the  Eu- 
clidean object  is  roughed  by  the  Gaussian  relationship).  The  geometry  of  the 
facets  and  the  direction  of  light  (assumed  to  be  from  an  infinitely  distant  source. 
so  all  rays  are  parallel^  determines  the  intensity  and  direction  of  specular  reflec- 
tion as  a  function  of  the  light  source  intensity,  the  normal  to  the  average  sur- 
face, the  direction  to  the  light  source  and  the  direction  to  the  viewpoint. 

This  model  has  to  be  modified  to  adapt  to  the  fractal  set  method.  In  the  fractal 
method,  there  is  no  need  to  rough  the  surface  to  provide  reflection  because  the 
surface  is  by  design  roughly  textured.  A  method  of  assigning  planar  front. 9  to 
each  pixel  space  has  to  be  determined  and  the  geometry  of  connecting  these 
fronts  identified.  With  these  modifications  to  the  lighting  model,  each  individual 
pixel's  color  intensity  can  be  modified  for  the  increase  in  intensity  associated  with 
the  light  which  falls  upon  it. 

The  model  also  allows  diffuse  reflection  (light  reflected  from  one 
object  to  another)  which  is  critical  to  bring  out  clarity  of  the  fractal  image.  For 
further  information  on  the  model  the  reader  is  referred  to  the  reference. 
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e.      Summary  of  the  Fractal  Mountain  Paradigm 

To  summarize  the  methodology  we  can  view  the  process  as  a  five 
step  process: 

-  Build  the  initiator  framework  or  stick  frame  model  of  the  mountain. 

-  Give  the  frame's  surface  texture  with  fractal  functions. 

-  Remove  hidden  surfaces  (pixels)  from  the  display. 

-  Illuminate  the  surface  with  lighting  algorithms. 

-  Project  the  surface  to  the  screen. 

B.    FRACTAL  TOOLS  FOR  TERRAIN  MODELING 

The  tools  presented  in  this  section  can  be  used  in  the  creation  of  fractal 
images  within  R3.  The  list  provides  a  basic  set  of  programming  tools  to  guide  the 
creation  process. 

1.     Equations  of  the  Lattice 

The  lattice  (or  controlling  structure)  can  be  very  useful  to  the  graphics 
programmer  to  implement  heuristics  or  bounding  functions  on  the  essentially 
random  progression  of  the  fractal  figure.  The  graphics  programmer  may  wish  to 
limit  the  growth  of  the  mountain  by  implementing  a  ground  level  plane  of  the 
lattice  and  a  maximum  height  that  the  mountain  can  obtain.  He  accomplishes 
this  by  arbitrarily  assigning  another  plane  of  the  lattice  as  the  upper  bounding 
plane.  The  height  of  the  mountain  can  then  be  checked  during  any  level  of  the 
fractal  recursive  descent  against  this  fixed  plane.  The  programmer  can  then  clip 
the  height  by  adjusting  the  random  equation  that  controls  the  upward  trend  to 
tend  towards  the  ground  again.  This  is  an  example  of  a  heuristic  applied  to  the 
fractal  recursion  that  controls  the  external  qualities  of  the  function. 

A  fractal  programmer  can  use  the  lattice  to  assign  the  initial  colors  to 
the  mountain  via  a  user  designed  set  of  rules.  The  lattice  aids  the  user  in  the 
implementation  of  the  rules  by  giving  reference  points  for  inclusion  of  branching 
conditions  (tree  line  to  snow  line  etc.)  and  can  be  used  in  conjunction  with 
decision  weights  to  add  a  varied  transition  between  textures.  The  determination 
of  the  initial  color  of  a  mapped  pixel  is  usually  a  controlled  random  process,  one 
of  the  primary  methods  of  control  being  the  lattice  or  some  derivative  thereof.  As 
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the  height  of  the  mountain  increases  (lattice  level),  it  becomes  increasingly  more 
likely  that  it  will  transition  to  another  texture.  This  can  be  controlled  by  adding 
the  lattice  level  as  a  factor  to  the  rule  that  decides  color. 

An  example  of  a  potential  lattice  equation  and  how  it  might  relate  to  the 
pixel  space  is  demonstrated  in  Figure  5.1.  The  actual  lattice  has  been  extended 
from  the  pixel  space  in  order  to  visually  demonstrate  how  it  relates  to  the  pixel 
set.  In  actuality,  this  is  not  the  case.  The  lattice  coincides  with  the  boundaries  of 
the  pixel  space.  Although  the  lattice  can  have  a  one-to-one  relationship  with  the 
pixel  space,  this  defeats  the  purpose  of  the  lattice  (macro  control).  By  grouping 
cubic  sets  of  pixels  into  a  well-formed  relationship,  we  can  better  implement 
heuristics  and  bounding  functions. 

As  a  lattice  example,  consider  a  pixel  space  that  is  created  by  abstracting 
the  real  world  coordinate  space  for  our  mountain  as  described  below.  We  desire  a 
real  world  space  to  be  a  cubic  area  established  by  the  box  20.000  ft. 
(x  coordinate  )  by  15,000  ft.  (y  )  by  20,000  ft.  (z  )  This  can  be  sectioned  into  a 
lattice  by  establishing  the  increment  of  distance  between  adjacent  lattice  points 
to  be  1000  ft.  and  establishing  the  corner  lattice  point  as  (0.0, 0)30. 

The  mapping  function  between  the  lattice  and  pixel  space  is  then 
straightforward.     The    size    of   the    pixel    (recall    equation    from    Chapter    3)    is 

: -=20  ft.     and    the     lattice     cubic     sections     contain     108  cubic  feet     or 

1,000  ft. 

equivalently  75,000  cubic  pixels. 

The  ground  level  can  then  be  identified  as  the  2000  foot  level  and  the 
bounding  height  can  be  assigned  a  level  of  10500  feet.  If  we  wish,  we  can  make 
the  bounding  heuristic  more  realistic  by  sectioning  the  lattice  into  mountainous 
areas,  each  having  different  bounding  levels. 


A  completely  arbitrary  set  of  dimensions,  increments  and  points. 
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The   Lattice    Control    Structure 


LATTICE 


FRACTAL  MOUNTAIN  STRUCTURE 


PIXEL  SET 


Fig  5.1.  The  Lattice  Control  Structure 
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2.  A  Fractal  Function  for  Contouring  Mountains 

The  usual  method  for  contouring  mountains  uses  a  randomized  variation 
of  the    mid-point    displacement    method    introduced    in    Chapter    4.    The    planar 

structure  is  typically  the  triangle31  imbedded  in  R3.  The  basic  methodology  is 
demonstrated  in  Figure  5.2.  Figure  5. 2. a  shows  a  triangle  with  its  first  iteration  of 
mid-point  displacement.  This  process  continues  until  all  triangles  have  reached 
the  desired  level  of  precision.  One  completed  structure  is  demonstrated  in  Figure 
5.2.b.  The  precision  is  typically  lower  (pixel  level)  than  that  demonstrated  in 
Figure  5.2.b  but  it  was  terminated  at  a  higher  level  to  better  demonstrate  the 
idea.  Random  techniques  (described  below)  are  used  to  produce  the  relatively 
accurate  picture  of  a  mountain  frame  as  depicted  in  Figure  5.2.c. 

In  practice,  the  random  techniques  are  implemented  with  the  mid-point 
displacement  function  during  the  fractal  recursive  descent.  The  random 
techniques  provide  local  disorder  to  the  fractal  function  which  provides  the 
computational  structure.  Results  have  shown  that  very  little  randomness  needs  to 
be  applied  to  the  regular  structure  of  Figure  5.2.b  to  achieve  satisfactory  results. 
The  mountains  created  for  the  film  Star  Trek:  The  Search  For  Spock  used  a 
limited  random  number  look-up  table  consisting  of  fewer  than  300  entries  [Ref. 

-]■ 

3.  The  Geometry  for  Mid-Point  Displacement 

The  general  approach  to  building  a  fractal  shape  as  illustrated  in  Figure 
5.2.c  is  to  use  the  algorithm  of  midpoint  triangle  displacement  combined  with  a 
randomized  displacement  along  the  normal  to  the  X-Z  plane  of  a  cartesian  three 
space  coordinate  system.  A  recursive  procedure  which  computes  this  relationship 
requires  as  inputs  the  points  of  the  triangle.  It  computes  the  midpoints  of  each 
line  of  the  triangle  and  inscribes  a  triangle  inside  of  the  initiating  triangle  by 
connecting  each  midpoint,  Figure  5. 3. a.  This  process  yields  four  triangles 
coincident  with  the  plane  of  the  initiating  triangle.  When  we  fix  the  X-Z  normal 
at   any  of  the  midpoints,  we  can  displace  the  midpoint   by  a  discrete  distance 


1    Any   regular  structure   suffices;   the   triangle   is  easy    to   use   and    yields   very    satisfactory 
results. 


71 


The    Triangular   Midpoint 
Displacement   Technique 


Determine 
Midpoints 


Connect 
Ml     to   M2 


Connect 
M2    to    M3 


Connect 
M3    to    Ml 


Fig  5. 2. a  The  1st  Iteration  of  Midpoint  Displacement. 


Fig  5.2.b  Completed  Structure. 


The  Random  Structure  was  rotated 
-30  degrees  around  the  X  axis  to 
accentuate    its    texture . 


Fig  5.2.c  Randomized  Version 


Figure  5.2.  The  Triangular  Midpoint  Displacement  Technique 
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along  the  normal  and  determine  a  point,  Figure  5.3.b.  Since  the  normal  is  to  the 
X-Z  plane,  it  is  sufficient  to  simply  modify  the  Y  coordinate  according  to  a 
positive  or  negative  value.  This  is  equivalent  to  displacing  the  midpoint  along 
the  X-Z  normal  up  or  down.  We  perform  this  displacement  to  each  midpoint 
normal  and  replace  the  midpoint  with  these  new  points.  This  yields  a  new 
structure  that  still  consists  of  four  triangles  but  with  each  coincident  with  a 
different  plane,  Figure  5.3.c. 

a.  Midpoint  of  a  Line  in  R3 

The  determination  of  the  midpoints  of  the  lines  of  the  initiating 
triangle  is  a  simple  process  that  uses  the  equation  of  Chapter  4.  and  fixes  the 
generator  ratio  constant  at  1.  This  simplifies  the  general  equation  of: 

Xj  +  (Generator. ratio. constant  x  X2) 
1  +  Generator. ratio. constant 
to  the  well-known  midpoint  relationships  of: 

X,  +  x2 

Xmid 
Y 

*  mid 

Zmid 

The  above  equations  completely  determine  the  midpoints  of  the  lines  formed  by 
each  endpoint  of  the  initiating  triangle. 

b.  Displacement  along  the  X— Z  Normal 

The  process  of  displacing  the  midpoint  along  the  X-Z  normal  is  a 
simple  one.  We  need  a  factor  such  that  the  displacement  can  obtain  a  varied 
magnitude.  This  is  best  aided  by  the  inclusion  of  a  random  variable  as  a  multiple 
of  some  scaling  factor  that  is  added  to  the  Y  coordinate  of  each  computed 
midpoint.  This  process  is  demonstrated  in  the  following  code  segment: 

Randvar     =  getrand(Seed); 

Point l[y]  =  Pointl[y]  +  (Scale  *  Randvar); 


2 

Y, 

+  Y2 

2 

z, 

+  z2 
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Four  Triangles  are 
created  for  each 
triangle  initiator 


Figure  5. 3. a.  Triangular  Midpoint  Displacement 


Computed    displacement 


After    1st    displacement 


X-Z  plane 


Figure  5.3.b  Displacement  Along  the  X-Z  Normal 


Negative    Y  displacement 


Original    generating    triangle 


Positive    Y  displacement 

Figure  5.3.c  Completed  Random  Fractal  Triangle. 
Figure  5.3.  The  Random  Midpoint  Displacement  Technique 
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A  valid  question  is,  why  the  normal  to  the  X-Z  plane?  There  are 
three  good  answers  to  this  question.  Using  the  normal  to  a  fixed  plane  simplifies 
the  computation  (eliminates  the  need  to  perform  planar  computations  at  each 
recursive  division).  It  also  is  generally  the  direction  that  we  want  the  mountain 
to  grow.  The  most  important  reason  however,  is  related  to  the  gapping  problem 
(described  below).  With  a  fixed  direction  for  displacement,  there  is  no  need  to 
communicate  the  direction  of  displacement  along  the  normal  between  adjacent 
side  computations.  The  recursive  levels  that  compute  adjacent  sides  are 
functionally  discordant.  It  is  demonstrated  below  that  the  solution  to  the  gapping 
problem  (inconsistent  random  numbers)  which  creates  the  need  to  communicate 
along  discordant  recursive  levels  is  algorithmically  difficult  to  solve  and  thus 
should  be  avoided. 

c.      The  Gapping  Problem 

One  problem  exists  for  the  midpoint  displacement  procedure  which 
utilizes  a  random  displacement  along  the  X-Z  normal  line.  It  is  indirectly  caused 
by  the  data  locality  aspect  of  the  inductive  process  of  the  recursive  fractal 
descent.  The  problem  exists  when  two  adjacent  sides  of  two  adjacent  triangles  are 
not  displaced  with  the  same  value.  Each  side  is  computed  during  independent 
levels  of  the  recursive  descent  so  there  is  no  practical  method  to  communicate  the 
random  numbers  for  the  displacement. 

The  gapping  problem  is  illustrated  in  Figure  5.4.  For  the  two 
triangles  that  are  extrapolated  from  the  structure,  there  is  an  unknown 
relationship  that  is  the  refndom  variable  used  to  displace  the  common  midpoint. 
If  triangle  A  uses  Rand  =  0.3  and  triangle  B  uses  Rand  =  -1.07,  then  the 
displacement  for  each  adjacent  midpoint  (which  are  at  the  start  coincident)  is 
skewed  in  the  opposite  direction.  This  creates  a  gap  in  the  fractal  landscape  that 
will  (in  all  likelyhood)  not  be  filled  by  other  fractal  shapes  from  neighboring 
triangles.  We  need  an  algorithm  which  can  insure  that  each  midpoint  (which  is 
always  shared  by  two  triangles)  has  the  same  displacement  along  the  normal  to 
the  plane. 
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d.     Solving  the  Gapping  Problem  via  Random  Tables 

The  solution  to  the  gapping  problem  is  straightforward  if  the 
programmer  adopts  the  random  number  table  as  his  random  function 
implementation.  The  goal  is  to  match  adjacent  triangles  with  a  seed  or 
displacement  within  the  random  table  so  that  the  random  number  returned  is 
equivalent  for  each  coincident  midpoint. 

There  exists  a  symmetry  within  the  triangle  of  Figure  5. 5. a  that 
allows  such  an  approach.  Ideally  we  want  the  point  Ma  to  be  displaced  by  the 
same  magnitude  when  triangles  T2  and  T4  (highlighted  by  textures)  compute 
their  random  numbers  for  Ma.  This  can  be  facilitated  by  the  inclusion  of  a  table 
seed  for  each  recursive  call  to  the  midpoint  displacement  routine  and  by  rotating 
the  orientation  of  the  midpoint  triangle  (the  triangle  created  by  the  three 
computed  midpoints)  labeled  T4  in  Figure  5. 5. a.  This  rotation  is  performed  in 
relation  to  the  random  table  and  not  in  relation  to  the  Cartesian  space.  It  is 
accomplished  by  adjusting  the  order  of  the  points  in  the  recursive  call. 

The  order  of  the  points  for  triangles  Tj  T2  and  T3  are  as  described 
in  Figure  5.5.b  and  for  T4  as  described  in  Figure  5.5.c.  All  four  triangles  generate 
a  recursive  sequence  and  use  the  same  seed  to  the  random  number  table.  The 
random  numbers  retrieved  from  the  table  must  observe  the  order  of  assignment 
that  is  demonstrated  in  Figure  5.5.d.  For  example,  the  line  segment  formed  by 
the  first  two  points  (Pj  and  P2)  input  to  the  midpoint  displacement  routine 
determine  the  midpoint  Rj.  This  midpoint  is  assigned  the  random  displacement 
from  the  table  corresponding  to  the  entry  seed.  The  next  midpoint  retrieves  the 
table  entry  corresponding  to  seed  +  1  and  so  on. 

If  this  technique  is  followed  the  sequence  of  random  numbers  will 
match-up  as  demonstrated  in  Figure  5.5.e.  The  recursive  calls  correspond  to  the 
code  segment  in  Figure  5.6. 
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The    Gapping   Proble 


m 


Ini  tinting    Triangle 


ADJACENT  TRIANGLES  SHARING 
A  COMMON  MIDPOINT 


Gene  rated    random    structure 


common    midpo i nt 


AFTER  MIDPOINT  DISPLACEMENT 


Gap    created    by    discontinuous 
random    displ acements    along 
the    Y   axi s  . 


Shared    common    midpoint 


Figure  5.4  The  Gapping  Problem 
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Solving    the    Gapping   Problem 


t3 


Shared 

mm 

/  Tl  > 

Ip^iik 

■a               \ 

ton    midpoint 


Figure  5 . 5 . a 


Figure  5 . 5 . b 


Given    that    the    input    to    the    random    table    is    seed=i 

Figure  5 . 5 . d 


i+l 


Figure  5 . 5 . e 


Figure  5.5.  A  Solution  to  the  Gapping  Problem 
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Main() 

{ 

LOAD  THE  INITIATING  TRIANGLE 
Seed  =  1; 

frac_triangle(P  j  ,P2,P3,Seed) 

} 

frac_triangle(Pj,P7,P3,Seed) 

{ 

DETERMINE  DISTANCE  BETWEEN  ENDPOINTS  OF  AN  INITIATOR 

If  (DIST  <  Pixel. length) 

{ 
Plotj>ixel(); 

return; 

} 

COMPUTE  THE  MIDPOINTS  (M1.M2.M3) 

ADJUST  THE  Y  COORDINATE  FOR  Mj  Using  Randtable(Seed) 
ADJUST  THE  Y  COORDINATE  FOR  M2  Using  Randtable(Seed+l) 
ADJUST  THE  Y  COORDINATE  FOR  M3  Using  Randtable(Seed  +  2) 

Seed  =  Seed  +  3; 

/*  Triangle  T,  */ 

frac_triangle(M1.P2,M2,Seed) 
/*  Triangle  T2  */ 

frac_triangle(M3.M2,P3,Seed) 
/*  Triangle  T3  */ 

frac_triangle(P  j  .M,  .M3,Seed) 
/*  Triangle  T4  '/ 

frac_triangle(M2,M3,M1,Seed) 
} 

Figure  5.6.  An  Algorithm  for  the  Midpoint  Displacement  Technique. 


79 


4.     Random  (Stocastic)  Fractals 

One  common  complaint  about  computer  graphics  images  and  animations 
is  the  artificial  perfection  of  the  displayed  shapes.  Our  mind  subconsciously  rebels 
against  the  order  that  is  displayed,  our  expectations  about  the  rough  reality  of 
nature  are  not  satisfied.  The  use  of  randomness  in  generating  fractal  images  is 
necessary  to  approximate  the  observed  disorder  of  nature.  An  example  is  the 
Koch  curve.  Although  it  resembles  a  snowflake,  it  lacks  the  realistic  look  that 
experience  trains  our  eyes  to  see.  In  a  mathematical  sense,  the  Koch  curve  is 
beautiful;  as  an  approximate  to  nature  it  lacks  appeal. 

To  approximate  the  rough  texture  of  nature,  we  are  forced  to  modify  the 
well-behaved  mathematical  relationship  of  the  initiator^generator  in  a  controlled 
manner  to  add  variety  to  our  computed  image.  This  modification  is  usually  by 
the  inclusion  of  a  random  variable  into  the  control  structure  within  the  fractal 
equation.  The  random  variable  must  exhibit  restraint.  It  cannot  be  allowed  to 
vary  wildly  without  structure. 

One    of    the    most    appealing    random    functions    which    provides    very 

satisfactory  results  in  fractal  images  is  the  normal  distribution32.     The  normal 
distribution  (as  opposed  to  a  uniform  distribution)   approximates  the  expected 
local  disorder  in  nature  (at  least  experimentally). 
a.     The  Normal  (Gaussian)  Distribution 

The  normal  distribution  is  used  throughout  the  natural  sciences  for 
many  applications.  It  was  first  derived  as  an  empirical  result  of  the  observed 
error  about  a  true  value  that  normally  occurs  when  measurements  are  taken  of  a 
natural  event.  The  symmetry  that  was  observed  from  error  measurement  and 
sampling  suggested  that  there  was  a  natural  order  to  such  observations.  These 
empirical  results  spurred  natural  scientists  and  mathematicians  to  try  to  fit  a 
curve  to  the  observed  graph  that  behaves  as  probability  requires  (i.e.  the  sum  of 
the     area     under     the     curve     equals     unity).     Many     of    the     early     scientists 


32  Often  referred  to  as  the  Gaussian  distribution,  the  normal  distribution  is  the  standard  bell 
curve  to  which  every  student  is  accustomed. 
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referred  to  the  normal  distribution  as  the  law  of  error  in  deference  to  its  roots  in 
experimental  natural  science. 

Many   functional   characterizations   of  the   normal   distribution   were 

developed33,  but  credit  is  usually  attributed  to  Carl  Frederic  Gauss  [Ref.  9:pp.  1- 
11]  who  formulated  a  least  squares  approach,  published  in  1809  in  Theoria  Motus 
Corporum  Coelestium.  The  form  of  the  normal  distribution  was  not  finalized 
until  the  early  20     century. 

We  take  our  definition  of  the  normal  distribution  from  [Ref.  9:pp  18]. 
refer  to  Figure  5.7. 

Definition: 

The  probability  density  function  of  a  normal  random  variable  X  is  given  by: 
ft  is  the  mean,  a  is  the  standard  deviation 
and  a2  is  the  variance. 


f(x:MV 


"Ik 


•exp 


2a2 


where    -ao<  x  <  oq  -oo<  fi  <oq  and  ct>0 


The    Normal    Distribution 


SCALE: 

horz    vert 


Figure  5  7.  The  Normal  Distribution 


ss  Many  mathematicians  can  lay  claim  to  founding  the  normal  distribution,  most  notably:  Pi- 
erre Simon  de  Laplace  and  Abraham  de  Moivre    Ref.  9:pp.  11 
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This  definition  is  the  general  case  of  the  normal  distribution.  We  are 
interested  in  the  behavior  of  the  function  and  need  a  practical  way  to  determine 
a  random  number  that  we  can  use  in  the  parameter  of  the  normal  to  the  plane  in 
the  geometrical  relationship  described  above.  To  facilitate  this,  we  simplify  the 
general  normal  distribution  to  the  well  known  standard  normal  distribution, 
illustrated  in  Figure  5.8.  The  standard  normal  distribution  function  is  the  special 
case  where  /i  =  0  and  a  =  1.  This  reduces  the  general  equation  to  the  simplified 
equation: 


f.M  = 


i 


'2* 


■exp 


x" 

T 


The  above  functions  describe  the  behavior  of  a  normal  random 
variable.  We  need  a  function  that  returns  values  from  that  function  which  will 
observe  the  period  of  the  normal  distribution.  This  means  we  need  a  string  of 
real  numbers  over  an  assigned  range  about  a  mean  that  will  observe  the 
frequency  of  the  normal  distribution. 


SCALE: 

Wz  a  5vert 


The    Standard 
^Normal     (Gaussian) 
Distribution 


/*=  0.0 

(7=  1.0 


N.  (x)  =  7±=  e  2 


>/iT 


■99.37%- 


Figure  5.8.  The  Standard  Normal  Distribution 
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b.     Standard  Computer  Random  Functions 

Some  computer  systems  provide  a  random  number  generating 
function  which  observes  the  normal  distribution.  If  this  is  provided,  then  it  can 
be  used  directly  (after  scaling)  as  a  parameter  to  displace  the  Y  coordinate  in  the 
geometry  of  the  normal  to  the  midpoint  displacement  as  described  above. 

Many  computer  systems  only  provide  a  random  number  generating 
function  which  is  uniformly  distributed  over  an  interval  of  integers.  This  is  a 
pseudo- random  number.  Such  a  function,  when  given  a  seed,  will  produce  a 
sequence  of  numbers  distributed  over  the  fixed  interval  defined  by  that  system. 
The  interval  is  typically  proportional  to  the  maximum  integer  defined  in  the 
compilers  of  the  system.  A  normal  distribution  routine  must  then  be  defined  that 
transforms  the  uniform  random  numbers  into  random  numbers  which  behave 
according  to  the  standard  normal  distribution  function. 

There  exist  transformation  functions  that  take  a  uniform  random 
variable  distributed  over  the  interval  [0,1]   into  an  approximate  normal  random 

variable  over  —  oo  <  x  <  oo  .  This  requires  the  uniform  random  variable  to  be 
mapped  into  the  interval  [0,1]  and  then  transformed  by  the  normal 
approximating  function. 

To  transform  a  uniform  random  variable  distributed  over  an  interval 
[0,max_int]  while  maintaining  the  distribution  density,  requires  the  following  step: 


UNFI0>1]  = 


UNF-0ma^n( 
max  int 


One  commonly  used  function  that  transforms  uniform  random 
variables  into  normal  random  variables  is  found  in  [Ref.  9:pp  49].  This  function 
uses  two  uniform  variables  from  [0,1],  denoted  UNF,  and  UNF2.  and  computes 
two  normal  random  variables,  denoted  NORMj  and  NORM,. 

NORM,  -  >/(-2logeUNF1  cos(2*rUNF2) 


NORMo  =  v/(-2logeUNF1  sin(2?rUNF2) 


S4  This  is  how  most  standard  system  provided  computer  subroutines  perform  the  operation 
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Appendix    B    contains    a    C    UNIX    routine    that    implements    an    algorithm    to 
compute  the  um'/orm[0,l]  -»  normal  [— oc,+oo]    transformation. 

A  programmer  must  be  very  careful  when  dealing  with  random 
number  generators  from  standard  system  subroutines.  These  routines  vary  widely 
and  can  provide  good  to  barely  adequate  results.  When  the  normal 
transformation  routine  is  written,  the  programmer  must  verify  experimentally 
that  his  function  adequately  models  the  normal  distribution.  This  process  is 
illustrated  by  Figure  5.9.  Appendix  B  also  contains  experimental  results  which 
verify  the  transformation. 

The  purist  may  not  accept  the  results  displayed  in  Figure  5.9  as  an 
accurate  transformation  (there  appears  to  be  a  skew  to  the  negative  direction). 
We  must  remind  ourselves  that  we  are  trying  to  approximate  the  roughness  of 
nature  and  minor  random  skewness  will  not  deter  us.  If  the  programmer  demands 
a  better  approximation,  it  is  a  simple  process  to  expand  the  sample  space  of  the 
test  and  build  a  table  with  exact  proportions  by  selective  deletion  of  skew 
density. 

c.      Random  Functions  versus  Table  Driven  Methods 

The  application  of  a  random  modifier  in  the  midpoint  displacement 
technique  can  be  achieved  via  two  methods. 

-  By  invoking  the  above  function  iteratively  as  a  variable. 

-  Or  by  a  variable  returned  from  a  table  lookup  operation  from  a  random 
table. 

The  choice  of  which  method  to  use  depends  on  the  programmer's  application  but 
each  has  its  ramifications. 

In  general,  the  table  lookup  operation  is  considerably  faster  than  the 
functional  method  but  must  by  its  definition  limit  the  amount  of  randomness  it 
contains.  The  major  issue  however  is  the  need  to  reproduce  a  figure  under  some 
requirement  for  fixed  terrain.  This  issue  was  the  driving  force  for  Loren  Carpenter 
from  Lucas  Film  in  determining  that  he  needed  to  use  a  table  driven  method  to 
produce  the  planet  images  for  the  film  Star  Trek:  The  Search  for  Spock  [Ref.  7]. 
He  had  to  be  able  to  fix  a  space  where  the  images  of  the  actors  could  be  imposed 


84 


Experimental    Results    for    the 

Computer    Generated    Normal 

Distribution 


Sample    Space: 

500    Random    Events 
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Transformation  Equations 
From  Uniform  Distribution 
Over  [0,1]  to  the  Normal 
Distribution  Over  -oo<x<oo 


N  2  =y    (-2lnU1  )     *  cos(2KU2) 


N  2=   /    (-2lnUJ  )     *    sin(2tU  ) 
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Experimental  Reiulti 
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Figure  5.9.  Experimental  Results  for  a  Computer 
generated  normal  distribution  over  [0,1] . 
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onto  the  fractal  images  and  could  not  allow  the  fixed  space  to  change  with  each 
frame  computed.  This  is  the  major  advantage  of  the  table  method.  By  retaining 
a  seed  to  a  table  of  random  numbers,  you  can  reproduce  the  sequence  of 
displacements  along  the  normal  during  the  fractal  recursive  descent. 

When  you  consider  the  existential  qualities  of  randomness  you  are 
confronted  with  basic  questions  about  determinism  and  order  in  the  universe.  It 
is  not  at  all  clear  which  rules  chance.  In  any  case,  we  can  deceive  perception 
with  a  relatively  small  table  of  random  numbers. 

The  question  of  how  much  randomness  is  enough  to  provide  for  a 
visually  appealing  texture  is  not  completely  clear.  In  [Ref.  7]  Smith  demonstrates 
a  variety  of  shapes  computed  with  the  same  algorithm  of  Figure  5.6  using 
random  number  tables  of  different  sizes.  He  demonstrated  that  as  few  as  five 
numbers  can  suffice  to  provide  enough  local  disorder  to  give  the  viewer  the 
acceptable  texture  of  a  mountain.  If  the  mountain  segments  are  viewed  at  the 
correct  perspective  and  scale,  this  perception  is  clearly  felt.  A  trained 
mathematician  would  find  the  five  element  mountain  statistically  unappealing 
however.  A  true  stochastic  construction  requires  a  continuous  random  function 
rather  than  a  discrete  table  method.  As  long  as  the  goal  of  our  computations  is 
merely  to  deceive  the  graphics  viewer,  it  suffices  to  use  the  random  number  table. 
The  table  must  be  large  enough  to  provide  for  an  appealing  textural  perception. 
A  complete  C  program  that  computes  a  triangular  mountain  segment  using  the 
random  displacement  midpoint  technique  is  contained  in  appendix  C. 
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VI.    SHORT  CUTS  TO  MOUNTAIN  SHAPES 

Since  the  fractal  mountain  computation  (the  full  approach  with  hidden 
surfaces  etc.)  is  so  costly  in  terms  of  resources,  it  is  important  for  us  to  consider 
shortcuts  that  can  lessen  this  burden.  This  is  best  realized  by  utilizing  the  hidden 
surface  and  curve  fitting  capabilities  that  are  provided  on  some  advanced 
graphics  systems. 

Our  goal  is  to  match  the  well  known  bicubic  surface  procedures  with  the 
structure  computed  by  the  simple  fractal  algorithms.  This  is  best  accomplished 
by    modifying    the    triangular    midpoint    displacement    technique    and    using    a 

rectangle35  as  the  basic  geometric  building  block.  Most  of  the  cubic  surface 
algorithms  use  the  rectangular  structure  as  their  basis,  so  it  is  easier  to  adapt 
them  to  our  fractal  structure. 

When  the  fractal  algorithm  of  Figure  6.2  has  its  computations  terminated 
before  reaching  the  level  of  pixel  size,  it  yields  a  connected  rectangle  structure 
like  the  one  shown  in  Figure  6.3.  This  structure  is  a  connected  Euclidean 
structure  that  can  be  used  as  a  base  on  which  other  algorithms  can  be  applied. 
Cubic  equations  can  fill  the  polygons  to  an  arbitrary  precision  and  standard 
hidden  surface  algorithms  can  eliminate  the  hidden  sides  of  the  computed 
surfaces  Simple  lighting  algorithms  can  be  applied  to  the  computed  surface  to 
achieve  a  realistic  lighting  effect36.  This  is  how  Voss  and  Carpenter  created  their 
fractal  surfaces  in  [Ref.  7]. 

A.    RECTANGULAR  MIDPOINT  TECHNIQUE 

Modifying  the  triangular  midpoint  algorithm  of  chapter  5  is  a  straightforward 
process  that  introduces  no  new  mathematics  or  difficulties.  It  consists  of  a 
procedure  to  split  the  midpoints  of  each  side  of  the  rectangle  and  a  procedure  to 
find  the  center  of  the  rectangle.    From  these  five  points,  we  construct  four  scaled 


55  We  actually  use  a  non-planar  four  sided  polygon.  We  refer  to  the  basic  structure  as  a  rec- 
tangle to  simplify  the  terminology. 

6  Couraud  shading  for  example. 
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rectangles,  as  demonstrated  in  Figure  6.1.  The  five  shared  midpoints  of  the 
generated  rectangles  are  then  displaced  along  the  normal  to  the  X-Z  plane 
according  to  a  random  Gaussian  value.  This  process  is  exactly  the  same  as  for  the 
triangular  algorithm  of  chapter  5.  The  gapping  problem  still  exists  and  this 
requires  an  algorithm  to  rotate  the  rectangle  relative  to  the  random  number  table 
and  the  starting  seed  to  insure  that  adjacent  midpoints  are  displaced  relative  to 
the  same  random  number.  The  basic  methodology  is  displayed  in  Figure  6.1,  the 
algorithm  is  contained  in  Figure  6.2  with  sample  results  in  Figure  6.3. 
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Four  Rectangles 
are  created  for 
each  rectangu- 
lar initiator. 


Computed 
Midpoints 


Computed    displacement 


Y   axis 


'X-Z  plane 


Negative    Y 
displacement 


Positive    Y   displacement  Original    generating    rectangle 


COMPLETED  FRACTAL  RECTANGLE 


Figure  6.1.  The  Rectangular  Midpoint  Displacement 
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Main() 

{ 
LOAD  THE  INITIATING  RECTANGLE 

Seed  =  1; 

frac_rectangle(P  j  ,P2,P3,P4,Seed) 

} 

frac^,ectangle(P1,P2,P3,P4,Seed) 

{ 

DETERMINE  DISTANCE  BETWEEN  ENDPOINTS  OF  AN  INITIATOR 

If  (DIST  <  Pixel.length) 

{ 
Plot  joint  (); 

return; 

} 

COMPUTE  THE  MIDPOINTS  (M1,M2,M3.M4,MC) 

ADJUST  THE  Y  COORDINATE  FOR  M,  Using  Randtable(Seed) 
ADJUST  THE  Y  COORDINATE  FOR  M2  Using  Randtable(Seed+l) 
ADJUST  THE  Y  COORDINATE  FOR  M3  Using  Randtable(Seed+2) 
ADJUST  THE  Y  COORDINATE  FOR  M4  Using  Randtable(Seed+3) 
ADJUST  THE  Y  COORDINATE  FOR  Mc  Using  Randtable(Seed+4) 

Seed  =  Seed  +  5; 

/*  Rectangle  R2  */ 

frac^;ectangle(P1,M1,Mc,M4,Seed) 
/*  Rectangle  R2  */ 

frac_rectangle(M2,Mc,M1,P2,Seed) 
I*  Rectangle  R3  */ 

frac_rectangle(M2,Mc,M3,P3,Seed) 
/*  Rectangle  R4  */ 

frac_rectangle(P4,M3,Mc.M4,Seed) 
} 


Figure  6.2.  An  Algorithm  for  the  (Rect.)  Midpoint  Displacement  Technique. 
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Rectangular    Mountain   Fractal 


Generating    rectangle 


The    recursive    termination    event    was    1/2    inch 


Figure  6.3.  An  Example  of  the  Rectangular  Mountain  Fractal 
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B.    PARAMETRIC  CUBIC  SURFACES 

A  complete  description  of  parametric  cubic  surfaces  is  too  involved  to  be 
described  in  this  study.  The  theoretical  basis  of  cubic  curves  is  not  directly 
applicable  to  fractal  geometry.  For  a  complete  description  refer  to  [Ref.  8:pp. 
514-536].  If  the  reader  is  already  familiar  with  cubic  curves  and  their  derivations, 
he  can  skip  by  the  section  on  cubic  curves  to  the  section  that  details  the 
application  of  cubic  surfaces.  For  any  reader  who  has  not  been  exposed  to  the 
derivations  of  parametric  equations  which  yields  cubic  curve  computational 
engines,  it  is  recommended  that  he  read  the  following  section  so  that  he  may  gain 
insight  into  the  mathematics  of  cubic  surfaces.  Detailed  knowledge  of  cubic 
curves  is  not  a  prerequisite  to  the  successful  use  of  cubic  surface  fitting  engines 
with  respect  to  fractal  surfaces.  It  is  helpful,  however,  to  understand  the 
underlying  mathematics  whenever  canned  equations  are  used. 

1.     Cubic  Curves 

The  general  method  of  cubic  curves  has  as  its  basis  that  any  continuous 
curve  in  R  can  be  expressed  in  parametric  form.  This  form  relates  the  points 
x,y,z  with  a  parameter  t  such  that  as  t  varies  within  some  range  of  values  the 
equations  solve  for  unique  points  on  the  curve.  Specifying  two  endpoints  and  two 
control  points  of  a  segment  of  the  curve  allows  us  to  define  certain  constraints  to 
be  applied  to  the  parametric  equations.  These  constraints  allow  us  to  manipulate 
the  parametric  form  of  the  equations  to  yield  a  simple  vector  product  definition 
of  that  segment.  Once  this  vector  product  is  established,  we  can  solve  for  points 
on  the  curve  by  picking  discrete  values  of  t  and  solving  for  x. y,z  in  turn.  This 
yields  a  discrete  approximation  of  the  curve  that  can  be  as  precise  as  needed, 
a.      Parametric  Cubic  Equations  of  a  Curve 

A  parametric  cubic  curve  is  one  for  which  the  points  in  R3  (x,y,z) 
are  each  represented  as  a  third-order  (cubic)  polynomial  of  some  ;>arameter  t. 
Because  we   deal   with   a  finite   segment   of  a  curve,   we   limit   the   range  of  the 


t  may  vary  between  0  and  1  for  example. 
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parameter  t  to  the  range,  0  ^  t    ^  1.  This  yields  the  equations: 

x(t)  =  axt3  +  bxt2  +  cxt  +  dx 
y(t)  =  ayt3  +  byt2  +  cyt  +  dy 
z(t)  =  azt3  +  b,t2  +  e,t  +  dz 
Each  equation  can  be  expressed  as  a  vector  product  as  x(t)  is  below: 


x(t)  =  [t3  t2  t  1] 


This  vector  product  separates  the  distinct  parameters  of  the  parametric  equation 
into  the  unknown  coefficients  of  x(t);  [ax'Dxcx^x]  an(^  ^e  parameter  t  that  we 
wish  to  manipulate.  Through  this  separation,  we  are  able  to  manipulate  them  as 
algebraic  entities.  If  you  multiply  the  vector  product  out.  you  find  that  the 
vector  product  is  equivalent  to  the  parametric  equation  that  precedes  it.  Denote 
this  product  as  x(t)  =  TCX  where 


T  =  ft3  t2  t  1 


and 


The  vector  T  is  the  same  for  x(t).  y(t)  and  z(t). 

We  now  establish  constraints  (as  a  set  of  control  points)  for  the 
equation  x(t)  evaluated  at  the  bounds  of  the  range  of  the  parameter  /.  (i.e.  t=0 
and  t  =  l).  We  consider  four  equations  of  x(t)  and  its  first  derivative  x'(t)  where 
these  boundary  conditions  yield  four  known  points. 

x(t)=TCx;  when  evaluated  at=0.  ->  x(0)  =[000  l]Cx 
x(t)  =  TCx;  when  evaluated  at=l.  -»  x(l)  =   [l  1  1  l]Cx 
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and  since  the  first  derivative  of  x(t)  is: 

x'(t)  =  [3t2  2t  1  0]CX  =  T'CX 

x'(t)  =  T'Cx;  when  evaluated  at=0,  ->  x'(0)  =  [0  0  1  0]CX 
x'(t)  =  T'Cx;  when  evaluated  at=l,  -►  x'(l)  =  [3  2  1  0]CX 

We  now  have  four  equations  that  can  be  grouped  into  a  vector  product: 


x(0) 

0    0    0    1 

x(l) 

1111 

x'(0) 

0    0    10 

x'(l) 

3    2    10 

We  recognize  that  x(0)  and  x(l)  are  the  endpoints  of  the  curve 
segment  and  x'(0)  and  x'(l)  are  components  of  the  tangent  vector  at  the 
endpoints  (y'(t)  andz'(t)  are  the  other  components).  With  this  knowledge  we 
are  able  to  solve  the  left  hand  side  of  the  equation  above.  These  points  (that  we 

call  Pj  through  P4)  are  the  control  points  that  we  establish  for  curve  fitting38. 
For  a  given  curve  segment  the  control  points  are  fixed.  We  rewrite  the  equations 
above  with  respect  to  these  known  control  points: 


0 

0 

0 

1 

1 

1 

1 

1 

0 

0 

1 

0 

3 

2 

1 

0 

Denote  this  equation  as: 


Gx  =  MCX 


The  matrix  Gx  is  often  referred  to  as  the  geometry  of  the  cubic  curve  and  M  as 
the  basis. 

This  equation  has  the  4  by  1  row  vector  Cx  as  the  only  unknown. 
The    elements    of   the    Cx    vector    are    the    parameters    (ax,bx.cx.dx)    from    the 


If  we  establish  ourselves  as  servers  then  these  four  points  are  the  user's  input  to  our  rou- 


tine. 
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parametric  equations.  We  can  solve  this  equation  for  these  parameters  and 
establish  the  parametric  equations  with  the  only  unknown  being  the  parameter  /. 
The  parameter  t  can  be  discretely  varied  over  its  range  of  0  ^  t  ^1.  providing  a 
set  of  points  on  the  curve.  It  is  through  these  constraints  that  the  control  points 
control  the  parametric  equations  and  produce  an  equation  that  can  produce  a 
discretely  sampleable  curve  segment  in  three  space.  Solving  the  equation  for  Cx  is 
straightforward: 

CX=M1GX 

Substituting  Cx  into  the  equation  for  x(t)  yields39: 

x(t)  -  TM    'Gx 

Similar  arguments  yield  the  equations  for  y(t)  and  z(t): 

y(t)  -  TM-'Gy 
z(t)  -  TM'G, 

The  matrix  M_1  is  constant  for  all  three  equations  and  is  usually 
denoted  by  the  type  of  surface  that  it  relates  to  Bezier  —  Mb,  Hermite  —  Mh 
etc.  It  is  through  the  control  points  and  their  interaction  with  the  constraints 
that  the  models  Bezier,  B-spline,  Cardinal  Spline,  Ferguson  (Hermite  or  Coon's) 
surface  etc.  modify  the  parametric  equations  and  provide  different  curve  fitting 
engines. 

For  each  model,  the  matrix  Mmodej  is  constant  throughout  all 
computations.  To  use  the  model  requires  the  determination  of  the  control  points 
(in  conjunction  with  how  they  relate  to  the  curve)  and  a  vector  multiplication 
engine.  Since  vector  pipeline  computations  are  ideally  suited  to  computers,  this 
method  becomes  a  fast  technology  for  curve  fitting  with  an  intuitive  appeal  for  a 
programmer. 


We  have  just  determined  the  Hermite  model  equation  for  x(t). 
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b.     An  Example:  Bezier  Cubic  Curves 

We  consider  the  model  called  Bezier  [Ref.  8:pp.  514-536].  The  Bezier 
model  defines  the  position  of  the  curve's  endpoints  and  uses  two  other  points  (not 
on  the  curve)  which  define  tangents  at  the  curve's  endpoints  (by  the  line  segment 
joining  the  tangent  points  to  the  endpoints). 

The  matrix  M  is  derived  by  setting  the  following  constraints  (see 
Figure  6.4).    One  endpoint  of  the  segment  is  located  at  P^ 


x(0)=P1 


The  other  endpoint  is  located  at  P4: 

x(l)=P4 

The  line  segment  from  Pj  to  P2  defines  a  tangent  at  P1  such  that  x'(0)  relates 
to  the  points  Pj,P2  as  below: 

x'(0)  =  SfPr-PJ 

And  similarly  for  the  tangent  at  P4  defined  by  P3,P4: 

x'(l)  =  3(P4-P3) 
Solving  for  Cx  in  terms  of  Mb  yields  the  cubic  Bezier  matrix  as: 


1 

3 
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1 

0 

0 

0 

G, 


Hence  the  equation  for  x(t)  is: 


x(t)  =   [ts  tz  t  1] 


1 

3 
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Pi 
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P. 
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Bezier    Curve 


Tangent,    defined    by 
P   P       line    segment. 


Tangent    defined    by 

P   P  „  line    segment. 
4     <1 


Figure  6.4.  An  Example  of  a  Bezier  Curve 


The  process  of  creating  a  Bezier  curve  given  the  above  parametric 
cubic  engine  is  a  simple  process  of  computing  discrete  points  on  the  curve  by 
substituting  values  along  the  range  of  t  and  fitting  the  curve  by  connecting  each 
point  with  a  line  segment.  This  provides  an  approximation  to  the  curve  that  can 
be  processed  at  an  arbitrary  precision  by  incrementing  St  with  smaller  and 
smaller  lengths. 

The  process  of  shaping  a  curve  is  accomplished  by  increasing  or 
decreasing  the  two  endpoint  tangents  formed  by  the  four  control  points.  It  can  be 
viewed  intuitively  by  thinking  about  each  tangent  as  a  force  which  pulls  the 
curve  in  the  direction  of  the  tangent  until  the  force  from  the  other  endpoint 
overcomes  the  original  at  the  midpoint.  The  two  endpoint  tangents  work  against 
one  another  proportional  to  the  distance  of  St  from  each  endpoint. 
2.     Bezier  Surfaces 

Extending  the  above  method  to  cubic  surface  sections  is  accomplished  by 
adding   a   new    parameter   s  that   we   vary   from  0^5    ^  1    as   we   did  with   the 
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parameter  t  in  cubic  curves.  The  connection  between  cubic  curves  and  surfaces 
can  be  made  by  fixing  one  parameter  and  varying  the  other  over  its  range.  This 
yields  a  cubic  curve.  The  equation  is  of  the  form  x(s,t)  and  is  written  as: 

x(s,t)  =  ans3t3  +  a12s3t2  +  a13s3t  +  a]4s3 

T   ^41'      t   ®42        "*"   ®43       '     ®44 

Written  in  the  algebraic  form: 

x(s,t)  =  SCxTe 

where  S  =  [s3,s2,s,l],  T  =  [t3,t2,t,l]  and  Tfc  is  the  transpose  of  the  matrix  T. 

The  complete  algebraic  manipulation  of  the  equation  to  arrive  at  the 
equation  below  is  similar  to  the  curve  process  as  described  in  the  previous 
section.  Its  details  are  covered  in  [Ref.  8:pp.  524-536].  The  equation  for  a  Bezier 
surface  patch  is: 

x(s,t)  =  SMbQ^^T* 

where  Mb  is  the  same  matrix  as  in  the  curve  equation,  Mb1  is  its  transpose  and 
Qx  is  the  x  component  of  sixteen  control  points  of  a  surface  patch.  Bezier 
surfaces  are  intuitive  in  their  appeal  and  serve  the  fractal  rectangular  mountain 
well.  To  apply  the  technique  to  the  mountain  of  Figure  6.3  requires  the 
application  of  a  routine  that  takes  the  non-planer  four  sided  shape  of  a  computed 
initiator  and  develops  a  connected  sixteen  point  figure  as  illustrated  in  Figure  6.5. 
The  inclusion  of  a  Bezier  subroutine  at  the  recursive  termination  event  after  this 
figtire  is  developed  matches  the  sixteen  point  figure  with  a  smooth  curve.  To 
achieve  edge  continuity  requires  that  adjacent  sides  have  the  same  four  points  in 
proper  juxtaposition  in  the  sixteen  point  matrix.  This  is  also  demonstrated  in 
Figure  6.5.  Bezier  surfaces  guarantee  such  continuity. 
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Bezier    Surface    Patches 


- — — , _ — , — _ — . — , — _ — __ — , — . — . — __ — , 


p^***-*  Bi  m 


Shared    edge    points 


Qx  for  R,  = 


Qfor  R,  = 


Figure  6.5.  Matching  Bezier  Surface  Patches  to  the  Fractal 
Rectangular  Mountain  Structure. 
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VII.    CONCLUSIONS 

A.    DIRECTIONS  FOR  FURTHER  STUDY 

Fractal  geometry  as  an  area  of  research  is  very  new.  Because  of  this,  there  is 
a  great  need  for  refinement  and  exploration.  What  is  known  needs  to  be  refined 
into  a  set  of  workable  techniques  with  reasonable,  simple  terminology  as  its  root. 
The  areas  that  are  unknown  need  to  be  explored  intrepidly.  With  this  goal  in 
mind,  the  following  paragraphs  quickly  review  some  areas  of  prospective  research. 
The  reader  is  invited  to  explore  their  potential. 

1.  Development  of  New  Fractal  Functional  Methods 

The  current  tools  of  fractal  functions  are  tentative  and  limited  in  their 
ability  to  yield  insight.  New  applications  of  the  recursive  initiator- generator 
paradigm  are  waiting  to  be  discovered.  This  area  of  research  is  especially  good  for 
the  graphics  programmer  since  the  graphics  medium  is  currently  the  best  method 
for  fractal  experimentation.  As  these  new  functions  are  developed,  they  can  be 
shared,  yielding  a  glossary  of  modeling  functions  that  can  be  molded  into  a 
cohesive  theory40.  Related  to  this  is  the  need  to  develop  a  functional  language 
(within  the  language  of  mathematics)  of  fractal  geometry  to  aid  in  the 
communication  of  ideas  and  in  the  eventual  coalescence  of  the  theory. 

2.  Fractal  Lighting  Model 

The  current  state  of  the  art  in  computer  graphics  lighting  models  lacks  a 
complete   model  for  the   pixel  set   paradigm  that  was   introduced   in  Chapter  5. 

There  are  a  great  many  practical  applications41  which  demonstrate  successful 
lighting  techniques  but  no  published  model  exists.  This  indicates  a  piecemeal 
undisciplined  adaptation  of  the  Euclidean  based  lighting  models.  Ray  tracing 
techniques  look  promising,  as  does  an  adaptation  of  the  Torrance- Sparrow 
lighting  model  that  was  discussed  in  chapter  5.  A  good  pixel  set  lighting  model 
would  open  the  avenue  of  complex  terrain  modeling  to  a  much  wider  audience. 


4    Nature's  fractal  map? 

As  evidenced  by  the  fractal  pictures  that  have  been  published. 
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3.  Fractal  Music 

In     [Ref.     7]     Voss    demonstrates    the    application    of    fractal    recursive 

techniques  to  —  noise  and  has  produced  interesting  if  not  pleasing  tonal  results. 

It  is  safe  to  surmise  that  sound  is  a  roughly  textured  physical  phenomenon  and 
that  it  may  be  possible  to  create  or  decipher  sound  using  a  fractal  model.  Such  a 
discovery  would  aid  science  in  the  area  of  (rapid  )  speech  recognition. 

4.  Fractal  Computer  Graphics  Architectures 

It  is  clear  from  our  discussion  that  new  computer  architectures  need  to 
be  developed  to  support  the  pixel  set  paradigm  and  the  computational  aspects  of 
fractal  functions.  Such  special  architectures  require  parallel  processing  capabilities 
coupled  with  vast  memory  resources.  A  real-time  fractal  terrain  image  generator 
is  one  such  architectural  possibility. 

5.  A  Better  Fractal  Definition 

Fractal  geometry  is  currently  attaining  a  wide  audience.  Because  of  that. 
it  is  time  that  trained  mathematicians  tackle  the  problems  associated  with  the 
imprecise  and  unworkable  current  definition  of  fractal  sets42.  That  definition  uses 
competing  definitions  of  dimension,  each  of  which  is  somewhat  difficult.  A  new 
definition  could  be  based  on  a  fractal  set's  functional  or  statistical  qualities.  Such 
a  definition  scheme  must  provide  tools  to  further  its  workability. 


42  Sadly,  there  has  been  little  attention  from  the  mathematical  community,  although  that  is 
changing.  It  is  with  great  timidity  that  ones  accepts  fractal  geometry  without  such  scrutin> 
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B.    CONCLUSIONS 

Fractal  geometry  is  an  old  idea  that  has  found  a  new  application  with  the 
advent  of  computer  imaging  techniques.  Its  acceptance,  has  spawned  a  great  deal 
of  research  and  has  provided  a  new  tool  to  observe  nature  through  a  different 
perspective.  We  must  be  careful  to  insure  that  our  findings  are  in  fact  valid.  We 
also  must  begin  the  coalescence  of  the  many  techniques  that  have  been  developed 
in  order  to  control  the  growth  of  this  concept  and  to  attain  true  scientific 
acceptance.    Without  this  acceptance  the  theory  will  be  criticized  (validly)  as  an 

imprecise  and  unproven  idea    .  This  would  be  an  unfortunate  occurrence  because 
of  the  potential  that  fractal  geometry  possesses. 

It  is  the  hope  of  the  author  that  this  work  has  illuminated  the  subject  of 
fractal  geometry  and  that  it  will  aid  others  in  their  research.  The  purpose  and 
essence  of  fractal  geometry  is  based  on  simple  concepts.  The  reader  must  not  be 
overawed  by  the  current  literature  and  should  retain  his  perspective  with  a  mild 
dose  of  skepticism.  He  must  not  be  blinded  by  skepticism  though  as  the  potential 
of  fractal  geometry  has  not  yet  been  realized.  In  the  final  analysis,  we  expect 
that  even  the  skeptical  reader  will  discover  the  mathematical  beauty  and 
applicative  power  that  fractal  geometry  possesses. 


This  of  course  is  the  current  state  of  affairs  with  fractal  geometry. 
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APPENDIX  A:  FRACTAL  COMPUTATION  IN  R2 

The  first  routine  is  the  main  routine  which  initializes  the  data  for  the  Koch  curve 
generator  and  initiates  the  recursive  process  on  each  side  of  the  initiator  triangle. 
The  second  routine  is  the  recursive  subroutine  which  performs  the  generator 
replacement  until  the  recursive  termination  event  is  reached.  The  termination 
event  is  defined  by  the  precision  of  the  desired  output  medium. 

KOCH.C 


/* 

This  is  the  main  program  which  controls  the  initialization  of 
the  koch  generator  parameters  and  initiates  recursive  operations 
on  each  side  of  the  initiator  triangle. 

*/ 

/*  Global  generator  and  initiator  data  */ 
int        Generator  joints; 

/*  The  number  of  points  in  the  GENERATOR    */ 
double  Gen_angle[lO]; 

/*  The  angle  formed  between  init  joint  1 
and  gen  joint  */ 

double  Gen_ratio[lO]; 

/*  The  between  init  joint  1  to  gen  joint  and 
gen  joint  to  init _point2  */ 

double  Tan_theda[lO]; 

/*  The  tangent  of  the  angle  formed  between 
init  joint  1  and  gen  joint  */ 

double  Curjoint[20][2]; 

/*  Vertices  of  initiator  structure  */ 
int        Object  joints  jmb; 

/+  The  number  of  vertices  of  the  initiating  structure  */ 

^include  <math.h>        /*  Standard  UNIX  include  file  for  math  library    */ 

f define x  0 
$  define y  1 
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/♦  BEGIN  MAIN  PROGRAM  7 

mainQ 

I 

/*  Local  variables  */ 
int  I; 

f*   Initialize  global  variables  "/ 

/*  Initial  points  ofthe  INITIATORS  for  demo  */ 

Cur   point [0][x]    -  4.0; 

Cur  point [0][y]  -  3.0  +  sqrt(3.0); 

( 'ur   point  [  1 J  [x]        5.0; 
Cur  point[l][y]  -  3.0; 

Cur  point  [2]  [x|    -  3.0; 
Cur  point  [2]  [y]   -  3.0; 

/*  Remember  to  close  the  side  of  the  triangle  4/ 

Cur  point[3][x]    -  4.0; 

Cur  point [3][y]    -  3.0  t  sqrt(3.0); 

Object    puts    iimb         3; 
Generator  points   =  3; 

/'  Angle  (in  radians  formed  between  init  point  1  and  gen  point 
for  demo)       '  / 

Gen   angle[l]    ■  0.0; 

Gen  angle[2]      0.4712388; 
Gen  angle[3]      0.0; 

I*   Ratio  of  distance  between  init    point  1  and  gen   point(i)  and 
distance  between  gen    point  (i)  and  init    point 2 

( icn  rat  iol  l]  0.5; 
Gen  ratio[2]  -  1.0; 
Gen  ratiof3       2.0: 


it)  4 


/<■■•■  ' 

r    ' 

"  i 

r  v 


GENERATE.C 


/* 

This  subroutine  computes  the  generator  from  a  given  set 
of  points  in  R2  that  define  a  line  segment  which  is  the 
initiator.  The  routine  is  recursive  and  terminates  at  a  predefined 
precision  that  is  input  to  the  subroutine. 

7 

/*  External  global  generator  data;  defined  in  main  subroutine  */ 
extern  int       Generator  joints; 

/*  The  number  of  points  in  the  GENERATOR    */ 
extern  double  Gen_angle[lO|; 

/*  The  angle  formed  between  init  joint  1 
and  gen  joint  */ 

extern  double  Gen_ratio[lO]; 

/*  The  between  init  joint  1  to  gen  joint  and 
gen  joint  to  init  joint2  */ 

extern  double  Tan_theda[lO]; 

/*  The  tangent  of  the  angle  formed  between 
init  joint  1  and  gen  joint  */ 

^include  <math.h>  /*  Standard  math  include  file  for  UNIX  lib  */ 


/*  BEGIN  RECURSIVE  PROCESS  */ 

generate(Xl,Yl,X2,Y2.precision) 
/*  Parameter  variables  */ 

double  XI, Yl,X2,Y2, precision; 


{ 


/*  Local  variables  */ 

long     N; 

double  Parray[3][2]; 

double  G_point[lO][2],DIST; 

double  Slope _init, Slope  perp, Slope  gen; 

double  X  jerp,Y  jerp,b  jerp,b_gen,TEMP: 

double  ten_thousand, one, zero, minus  jne; 

int       I, J; 
/*  assign  constants  */ 

ten  thousand  =  10000.0;  one  =  1.0;  zero  =  0.0;  minus  one  =  -1.0: 
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/*  The  Koch  curve  is  defined  in  the  infinite  but  our  re<- 
will  terminate  after  the  distance  between  points  becomes 
then  the  length  of  a  pixel.    */ 

/*  Determine  distance  between  point  1  and  point  2  "/ 

TEMP  =  (X2  -  X1)*(X2  -  XI j  -    Y2  -  Yl  "  Y2  -  Yl 
DIST  =  sqrt(  TEMP    ; 

/*  IF  DIST  less  than  the  precision  then  terminate  this 
recursion  and  begin  backtracking     */ 

if  (DIST  <  precision) 

{ 

/*    Put  your  Point  plotting  routine  here  " 

printf( "polyline  2"): 
printf("cxf  %l  0.000000". Xl.Yl  : 
printf("c7.f  cwcf  0.000000". X2.Y2  : 
return; 

} 

/*  Put  INITIATOR  points  one  and  two  into  the  first  ar.  : 
points  of  the  GENERATOR  points  array  as  they  are  alw      - 
part  of  the  generated  structure  */ 

Gjoint'l[l]  =  XI; 

G_point[l][2J  =  Yl; 

G_point[Generator_points  +-  2j  lj  =  X2: 
G_point[Generator_points  -2   2=  Y2: 
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/•*  Determine  the  slope  of  the  line  formed  by  the  init _pointl 
and  init  joint2.  This  is  the  slope  of  the  INITIATOR  */ 

if  (X2  !=  XI) 

{ 

if  (Y2  !=  Yl) 

{ 
Slope Jnit  =  (Y2  -  Yl)/(X2  -  Xl); 

} 
else 

{ 

Slope_init  =  0.0; 

} 
} 

else 

{ 

/*  We  can't  have  infinity  in  a  register 

so  settle  with  10k  */ 

Slope_init  =  ten_thousand; 

} 

/*  For  each  GENERATOR  point  (except  end  points  as  they  are  equal 
to  the  INITIATOR  end  points)  find  the  X,Y  values.  This  is 
accomplished  by  using  the  data  from  the  global  external  variables. 
The  constant  data  about  the  ratios  and  angles  between  the 
INITIATOR  and  GENERATOR  remain  the  same  regardless  of  the 
INITIATORS  length  or  position  in  EUCLIDIAN  space  */ 

for  (1=1;  I  <=  Generator_points;  I+  +  ) 
{ 

/*  Using  the  ratios  of  the  generator  perpendicular  intercept 
points  on  the  INITIATOR  determine  the  X,Y  values  of  the 
point  of  intersection  of  the  perpendicular  from  the 
GENERATOR  point  to  the  INITIATOR  line.  */ 

Xjerp  =  (XI  +  Gen_ratio[I]  +  X2)/(1.0  +  Gen_ratio[I]); 
Yjerp  =  (Yl  +  Gen_ratio[I]  *  Y2)/(1.0  +  Gen_ratio[lj); 
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/*  If  the  angle  of  the  INITIATOR  point  1  and  the  GENERATOR 
point  in  question  is  zero  then  the  GENERATOR  point  is 
coincident  with  the  INITIATOR  line  and  no  further 
calculations  are  necessary  '/ 

if  (  Gen_angle[I]  ==  zero  ) 
G_point[I+l][l]  =  X_perp: 
G_point[I+l][2]  -  Yjperp: 

} 
else 

{ 

/*  There  are  three  STATES  possible  at  this  time.  STATE  1 
where  the  slope  of  the  initiator  line  is  parallel 
to  the  X  or  Y  axis  (which  causes  havoc  with  the  line 
equations).  STATE  2  where  the  slope  of  the  line  formed 
by.  the  initiator  point  1  and  the  unknown  generator  point 
is  parallel  to  the  X  or  Y  axis.  Or  STATE  3  where  no  lines 
are  parallel  to  any  axis.    */ 

/*  Determine  the  slope  of  the  line  through  the  INITIATOR 
point  1  and  the  unknown  GENERATOR  point  using  the 
tangent  of  the  Gen_angle  in  Init.h  */ 

Slope_gen    =  (Tan_theda[I]  +  Slope_init)/ 
(one  -  Tan_theda[I]  *  Slope_init); 

if  ((Slope_gen  !=  zero)  && 

(Slope_gen  <  ten_thousand  )) 

{ 

/*  Condition  one  of  STATE  3  */ 

/*  Determine  Y-intercept  for  the  generator  line  */ 

b_gen  =  Yl  -  (Slope_gen  *  Xl); 

if  ((Slope_init  ==  zero)  | 
(Slope_init  ==  ten_thousand)) 

{ 

/*  STATE  1  */ 
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if  (Slope_init  ==  ten_thousand  ) 

{ 
/*  STATE  1  condition  1;  INITIATOR  is  parallel 

to  the  Y  axis  */ 

G_point[I+l][2]  =  Yj>erp; 
G_point[I+l][l]  =  (G_point[I+l][2]  -  b_gen)/ 
Slope_gen; 

} 

else 

{ 
/*  STATE  1  condition  2;  INITIATOR  is  parallel 

to  the  X  axis  */ 

G_point[I+l][l]  =  X  j>erp; 
G_point[I+l][2]  =  Slope_gen  * 

G_point[I+l][l]  +  b_gen; 

} 

}  /*  END  STATE  1  */ 

else 

{ 

/*  STATE  3  */ 

/*  Determine  slope  of  perpendicular  line  through  the 
INITIATOR  perpendicular  intercept.  */ 

Slope_perp  =  (minus_one)/Slope_init; 

/*  Determine  Y-intercept  for  perpendicular  line  */ 

b_perp        =  Y_perp  -  (Slope  jerp  *  X_perp); 

/*  Determine  the  X,Y  values  of  the  unknown  GENERATOR 
point.  */ 

G_point[I+l][l]  =  (bjperp      -  b_gen        )/ 

(Slope_j;en  -  Slope  _j>erp); 
G_point[I+l][2]  =  Slope_gen  * 

G_point[I+l][l]  +  b_gen; 

} 
}  /*  END  STATE  3  cond.  1  if  */ 
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else 

{ 

/*  STATE  2  */ 

Slope_perp  =  (minus_one)/Slope_init; 
b_perp        =  Y^perp  -  (Slope j>erp  *  X_perp): 

if  (Slope_gen  ==  one) 

{ 
G_point[I+l][l]  =  XI; 

G_point[I+l][2]  =  Slope^perp  *  G_point[I+l][l] 
+  b_perp; 

} 
else 

{ 

G_point[I+l][2]  =  Yl; 

G_point[I+l][l]  =  (G_point[I+l][2|  -  b_perp)/ 
Slope_perp; 

} 

} 
}    /*  END  IF  7 

}    /*  END  FOR  7 

/*  Start  recursion  on  each  line  formed  by  the  generator  from 

right  to  left  */ 

for  (J=l;  J  <=  Generator_points  +  1;  J  +  +  ) 

{ 
generate(Gj?oint[J]    [l],G_point[Jj    [2], 

G_point[J  +  l][l],G_point[J+l][2].  precision): 

} 
/*  END  generate  */ 

} 
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APPENDIX  B:  RANDOM  NUMBER  GENERATORS 

The      routine      below      is      a      C      UNIX      UCB      implementation     of     the 
uniform  distribution   [0,1]  -*  standard    normal    [— oo,oo]  transformation.  It 

generates  a  500  entry  table  of  random  numbers  that  observes  the  period  of  the 
standard  normal  distribution.  Following  this  routine  are  statistics  that  verify  the 
transformation. 

RANDOM. TABLE. GENERATOR.C 


I* 

This  subroutine  will  build  a  table  in  memory  that  contains  500  random 
numbers  that  observe  the  period  of  a  standard  normal  variable 

7 


^include  <math.h>        /*  Standard  UNIX  include  file  for  math  library    */ 

/*  External  global  variables  */ 
extern  double  RAND [500]; 

/*  BEGIN  MAIN  PROGRAM  */ 

rand_table_gen  ( ) 

{ 
/*  Local  Variables  */ 

int       I,J; 

double  UNF1,  UNF2; 
double  range, pi; 
int       factor; 

pi  =  3.1415926535; 
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/*  Determine  the  range  for  the  random  numbers  of  UNIX  UCB  "/ 

range  =  2; 

for  (J=l;  J<=30;  J  +  +  ) 

{ 

range  =  range  *  2; 

} 
range  =  range  -  1; 

/*  Set  the  random  number  generator  seed  */ 

srandom(475836); 

/*  Create  a  Table  for  500  entries  */ 

for  (1=0;  I<  500;  1  =  1  +  2) 

{ 
/*  Get  a  uniform  random  number  through  the  Unix  C  subroutine  */ 

UNF1  =  randomQ; 
UNF2  =  randomQ; 

/*  Normalize  the  uniform  random  number  to  the  interval  [0,1]  */ 

UNF1  =  UNF1  /  range; 
UNF2  =  UNF2  /  range; 

/*  Mold  the  uniform  random  variable  into  the  approximate  normal 
distribution  */ 

factor  =  1.0; 

if  (log(UNFl)  <  0.0)  factor  =  -1.0; 

RAND[I]  =  sqrt(factor  +  (2.0  *  log(UNFl)))  * 

cos  ((2*pi*UNF2)); 
RAND[I]  =  RAND[I]  *  factor; 
factor  =  1.0: 

if  (log(UNF2)  <  0.0)  factor  =  -1.0; 
RAND[I+1]  =  sqrt(factor  *  (2.0  *  log(UXFl)))  * 

sin  ((2*pi*UNF2)); 
RAXD[I+1]  =  RAND[I+1]  *  factor; 

} 
return; 
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VERIFYING  STATISTICS 

The  UNIX  UCB  operating  system's  uniform  distribution  random  number 
generating  function  spans  the  interval  defined  by  its  integer  range.  For  a  VAX 
11/780  implementation  this  is  equivalent  to  231  -  1  or  [0,2147483647]. 

The  random  number  seed  was  assigned  the  value  of  475836.  The  UNIX  UCB 
random  number  generator  with  a  fixed  seed  yields  a  fixed  sequence  of  numbers 
returned  from  the  function,  uniformly  distributed  over  the  range.  This  yields  a 
valuable  function  if  the  table  needs  to  be  reproduced  with  the  same  sequence 
after  transformation. 

The  table  below  shows  the  results  of  the  uniform  distribution  sequence  after 
it  was  squeezed  into  the  interval  [0,1].  These  results  show  that  the  uniform 
distribution  has  an  acceptable  distribution  over  its  range.  The  transformation 
into  [0,1]  preserves  the  distribution  from  the  original  range  ([0,231  -  l]). 

Analysis  of  the  normalized  uniform  random  numbers 

0.0  -*  0.1  =  52 

0.1  ->  0.2  =  47 

0.2  -►  0.3  =  44 

0.3  ->  0.4  =  49 

0.4  -+  0.5  =  49 

0.5  -*  0.6  =  47 

0.6  ->  0.7  =  51 

0.7  -*  0.8  =  57 

0.8  -  0.9  =  48 

0.9  ->  1.0  =  56 

The  table  below  shows  the  distribution  after  the 

uniform  distribution  [0,1]  -*  standard  normal  [— oo,oc]  transformation  given  the 
numbers  as  described  in  the  above  table.  This  is  the  data  which  was  used  to 
build  Figure  5.9.  The  transformation  is  acceptable  for  the  purpose  intended,  that 
is,  to  simulate  nature's  perceived  disorder  in  a  fractal  function. 
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Analysis  of  the  normal  (Gaussian)  random  numbers 

X  <=  -2.75  =  5 

-2.75  <  X  <=  -2.25  =  8 

-2.25  <  X  <=  -1.75  =  16 

-1.75  <  X  <=  -1.25  =  29 

-1.25  <  X  <=  -0.75  =  56 

-0.75  <  X  <=  -0.25  =  88 

-0.25  <  X  <=    0.25  =  115 

0.25  <  X  <=    0.75  =  87 

0.75  <  X  <=    1.25  =  59 

1.25  <  X  <=    1.75  =  21 

1.75  <  X  <=    2.25  =  12 

2.25  <  X  <=    2.75  =  4 

2.75  <  X  =0 


115 


APPENDIX  C:  THE  TRIANGULAR  MOUNTAIN 

The  first  routine  is  the  main  routine  which  initializes  the  generator  data  for 
the  initiating  triangle  and  initiates  the  recursive  process.  The  second  routine  is 
the  recursive  subroutine  which  performs  the  generator  replacement  until  the 
recursive  termination  event  is  reached,  which  is  defined  by  the  precision 
parameter. 

MOUNTAIN. C 


/* 

This  is  the  main  program  that  controls  the  initialization  of  the  triangle 
initiating  structure  and  initiates  the  recursion  on  that  triangle.  The 
recursion  will  proceed  until  the  recursive  termination  event  (defined 
by  the  precision  global  parameter) 

*/ 


^include  <math.h>        /*  Standard  UNIX  include  file  for  math  library    */ 

/*  Global  Defines  */ 

#  define  x  0 
f  define  y  1 
^define  z  2 

/*  Global  Tables  */ 

double  RAND[500]; 
double  Precision; 
double  Scale; 
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/*  BEGIN  MAIN  PROGRAM  */ 

mainQ 

{ 

/*  Local  Variables  */ 

int       I,J,K; 

double  P1[3],P2[3],P3[3]; 

int       Seed; 

/*  Create  the  initiating  triangle  structure  */ 

Pl[x]  =  4.5; 
Pl[y]  -  3.25; 
Pljz]  =  0.0; 

P2[x]  =  7.0; 
P2[y]  =  3.25; 
P2[z]  -  0.0; 

P3[x]  =  5.75; 

P3[y]  =  3.25  +  sqrt(((2.5  *  2.5)  -  (1.25  +  1.25))); 

P3[z]  =  0.0; 

/*  Build  the  random  number  table  (appendix  B)  */ 

rand_table_gen(); 

/*  Fractalize  until  desired  precision  */ 

Seed  =  0;  /*  Entry  seed  to  the  random  number  table  */ 

Precision  =  0.3;    /*  Recursive  termination  distance  */ 
Scale  =  0.2;  /*  Scaling  factor  for  vertical  Y  displacement 

in  the  mountain_generate  subroutine  */ 

mountain_generate(Pl,P2.P3.Seed); 

/*  END  MAIN  */ 
} 
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GENERATE  MOUNTAIN. C 


/* 

This  is  the  subroutine  that  computes  the  four  generated  triangles  from 
an  initiating  triangle.  The  routine  is  recursive  and  terminates  at  a 
predefined  precision  defined  in  the  global  parameter  Precision 

r 

^include  <math.h>  /*  Standard  math  include  file  for  UNIX  lib  */ 

f  define  x  0 
#■  define  y  1 
#  define  z  2 

/*  Global  Structures  */ 

extern  double  RAND[500]; 
extern  double  Precision; 
extern  double  Scale; 

/*  BEGIN  RECURSIVE  PROCESS  */ 

mount  ain_generate(Pl,P2,P3, Seed) 

/*  Parameter  variables  */ 

double  P1[3],P2[3],P3[4]; 
int       Seed; 


/*  Local  variables  */ 

int       I,J; 

double  Pgenl  [3] ,Pgen2 [3]  ,Pgen3 [3] ; 
double  Pmidl[3],Pmid2[3].Pmid3[3]; 
double  TEMP,DIST,TWO; 

TWO  =  2.0; 
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/*  Determine  distance  between  point  1  and  point  2  * / 

TEMP  =  (P2[xj  -  Pl[x])  +  (P2[x]  -  Pl[x|)  + 
(P2[y|-Pl[y])*(P2[y]-Pl[y])  + 
(P2[z]  -  Pl[z])*(P2[z]  -  Pl[z]); 
DIST  =  sqrt(  TEMP  ); 

/*  IF  DIST  less  than  one  then  terminate  this  recursion  and 
begin  backtracking  */ 

if  (DIST  <  Precision) 

{ 
/*  Put  your  polygon  plotting  routine  here  */ 

printf  ( "poly  gon3  " ) ; 

printf(MO/of  %f  %f\Pl[x],Pl[y],Pl[z]): 

printf("%f  %f  %f \P2[x],P2[y],P2[zj): 

printf("%f  %f  %P',P3[x],P3[y],P3[z]); 

return; 

} 
/*  Manage  the  Seed  number  for  a  500  entry  table  */ 

if  (Seed  >  496)  Seed  =  0; 

/*  Find  the  midpoints  of  each  triangle  leg  */ 

for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x.y,z  */ 

{ 

Pmidl[I]  =  (P1[I]  +  P2[I])  /  TWO; 

} 

for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x,y,z  */ 

{ 
Pmid2[I]  =  (P2[I]  +  P3[I|)  /  TWO; 

} 

for  (1=0;  I<  =  2:  I+  +  )    /*  0  thru  2  =>  x.v.z  */ 

{ 

Pmid3[I]  =  (P3[I]  +  P1[I])  /  TWO; 

} 
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/*  Adjust  the  Y  coordinate  =>  normal  from  Z-X  plane  */ 

Pmidl[y]  =  (Scale  *  RAND[Seed])      +  Pmidl[y]; 
Pmid2[y]  =  (Scale  *  RAND[Seed+l])  +  Pmid2[y]; 
Pmid3[y]  =  (Scale  *  RAND[Seed+2j)  +  Pmid3[yj; 
Seed        =  Seed  4-  1; 
/*  Recurse  on  the  triangles  according  to  the  reverse  order  rule 
for  the  interior  triangle  to  preserve  seed  order  */ 

mountain_generate(Pmidl,P2,     Pmid2,Seed); 
mountain_generate(Pmid3,Pmid2,P3.     Seed); 
mountain_generate(Pl,     Pmidl,Pmid3,Seed); 
mountain_generate(Pmid2,Pmid3,Pmidl,Seed); 
/*  END  generate  */ 

} 
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APPENDIX  D:  THE  RECTANGULAR  MOUNTAIN 

The  first  routine  is  the  main  routine  which  initializes  the  generator  data  for 
the  initiating  rectangular  shape  and  initiates  the  recursive  process.  The  second 
routine  is  the  recursive  subroutine  which  performs  the  generator  replacement 
until  the  recursive  termination  event  is  reached,  which  is  defined  by  the  precision 
parameter. 

RECTANGULAR-MOUNTAIN. C 

/* 

This  is  the  main  program  that  controls  the  initialization  of  the  rectangular 
initiating  structure  and  initiates  the  recursion  on  that  rectangle.  The 
recursion  will  proceed  until  the  recursive  termination  event  (defined 
by  the  precision  global  parameter) 

7 

^include  <math.h>        /*  Standard  UNIX  include  file  for  math  library    */ 
/*  Global  Defines  */ 

fdefine  x  0 

#  define  y  1 

#  define  z  2 

/*  Global  Tables  */ 

double  RAND [500]; 
double  Precision; 
double  Scale; 
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/*  BEGIN  MAIN  PROGRAM  */ 

main() 

{ 
/*  Local  Variables  */ 

int      I,J,K; 

double  P1[3],P2[3],P3[3],P4[3]; 

int       Seed; 

/*  Create  the  four  sided  polygon  initiating  structure  */ 
Pl[x]  =  2.0; 
Pl[y]  =  5.0; 
Pljz]  =  1.0; 

P2[x]  =  6.0; 
P2[y]  =  5.0; 
P2[z]  =  1.0; 

P3[x]  =  7.5; 
P3[y]  =  6.5; 
P3[z]  =  3.0; 

P4[x]  =  3.5; 
P4[y]  =  6.5; 
P4[z    -  3.0; 

/*  Build  the  random  number  table  (appendix  B)  */ 
rand_table_gen(); 

/*  Fractalize  until  desired  precision  */ 
Seed  =  100; 
Precision  =  0.5; 
Scale  =  0.07; 

mountain_generate(Pl,P2,P3,P4.Seed); 

/*  END  MAIN  */ 
} 
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RECTANGULAR. GENERATE. MOUNTAIN. C 


/* 

This  is  the  subroutine  that  computes  the  four  generated  rectangles  from 
an  initiating  rectangle.  The  routine  is  recursive  and  terminates  at  a 
predefined  precision  defined  in  the  global  parameter  Precision 

l* 

^include  <math.h>  /*  Standard  math  include  file  for  UNIX  lib  */ 

§  define x  0 
$  define  y  1 
#  define  z  2 

/*  Global  Structures  */ 

extern  double  RAND [500]; 
extern  double  Precision: 
extern  double  Scale; 

/*  BEGIN  RECURSIVE  PROCESS  */ 

mountain_generate(Pl,P2,P3,P4,Seed) 

/*  Parameter  variables  */ 

double  P1[3],P2[3].P3[3].P4[3); 
int       Seed: 


{ 


/*  Local  variables  */ 
int       I,J; 

double  Pmidl[3].Pmid2[3],Pmid3[3].Pmid4[3],Center[3]; 
double  TEMP,DIST.T\VO.FOUR: 

TWO  =  2.0:  FOUR  -  4.0: 

/*  Determine  distance  between  point  1  and  point  2  */ 
TEMP  =  (P2[x]  -  Pl[x)nP2[x]  -  Pl[x]j  + 
(P2[y]  -  Pl[y]r(P2[y]  -  Pl(y])  + 
(P2[z]  -  Pl[z])*(P2[z]  -  Pl[z]); 
DIST  =  sqrt(  TEMP  ); 
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/*  If  DIST  less  than  one  then  terminate  this  recursion  and 
begin  backtracking  */ 

if  (DIST  <  Precision) 

{ 
/*  Put  your  Polygon  output  routine  here  */ 

printf  ( "polygon4  " ) ; 

printf("%f  %f  %f",Pl[x],Pl[y],Pl[z]); 

printf("%f  %f  %f\P2[x],P2[y],P2[zj): 

printf("%f  %f  %f\P3[x],P3[y],P3[zj): 

printf("%f  %f  %f*,P4[x],P4[y],P4[z]): 

return; 

} 

/*  Manage  the  Seed  number  for  a  500  entry  table  */ 
if  (Seed  >  496)  Seed  =  0; 

/*  Find  the  midpoints  of  each  rectangle  leg  */ 
for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x,y,z  */ 

{ 
Pmidl[I]  =  (P1[I]  +  P2[I])  /  TWO; 

} 

for  (1=0;  I<  =  2;  I+  +  )    /*  0  thru  2  =>  x,y,z  */ 

{ 

Pmid2[I]  =  (P2[I]  +  P3[I])  /  TWO; 

} 

for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x,y,z  */ 

{ 
Pmid3[I]  =  (P3[I]  +  P4[I])  /  TWO; 

} 

for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x.y,z  */ 

{ 
Pmid4[I]  =  (P1[I]  +  P4[I])  /  TWO: 

} 
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/*  The  four  sided  polygon  is  non-planar  so  average  the  xyz-displacement 
for  a  best  fit  approach  */ 
for  (1=0;  I<=2;  I+  +  )    /*  0  thru  2  =>  x.v.z  '/ 

{ 
Center[I]  =  (P3[I]  +  PI  [I]  +  P4[I]  +  P2[I])  /  FOUR; 

} 

/*  Adjust  the  Y  coordinate  =>  normal  from  Z-X  plane  */ 
Pmidl[y]  =  (Scale  *  RANDjSeed])  +  Pmidl[y); 
Pmid2[y]  =  (Scale  *  RAND[Seed+l])  +  Pmid2[y]; 
Pmid3[y]  =  (Scale  *  RAND[Seed+2])  +  Pmid3(yj; 
Pmid4[y]  =  (Scale  *  RAND[Seed+3])  +  Pmid4[yj; 
Center[y]  =  (Scale  *  RAND[Seed+4])  +  Center[y]; 
Seed         =  Seed  +  4; 

/*  Recurse  on  the  rectangles  according  to  the  reverse  order  rule 
for  the  interior  rectangles  to  preserve  seed  order  */ 

mountain_generate(Pl,  Pmidl.  Center, Pmid4, Seed) 
mountain_generate(Pmid2, Center, Pmidl,  P2,  Seed) 
mountain_generate(Pmid2, Center, Pmid3,  P3,  Seed) 
mountain_generate(P4,     Pmid3,  Center, Pmid4, Seed) 

/*  END  generate  */ 
} 
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APPENDIX  E:  GEOMETRIC  SUPPORT 

Many  fractal  applications  and  computer  graphics  models  use  the  normal  to  a 
plane  as  a  computational  reference  point.  For  this  reason,  this  appendix  is 
devoted  to  two  tools  for  determining  the  plane  equation  of  a  polygon  and  the 
equation  of  the  normal  to  the  computed  plane. 

Determinant  Approach  to  the  Planar  Equation 

One  of  the  most  common  forms  of  a  planar  equation  is  the  general  form.  This 
form  uniquely  describes  a  plane  through  four  coefficients  A,B,C  and  D: 


Ax   +  By   +  Cz   =  D 

With  three  points  on  a  plane,  you  can  determine  the  planar  equation  by 
computing  the  coefficients.  This  approach  utilizes  the  determinant  form  of  the 
planar  equation.  Given  the  points  Pj  =  (xl,y  1,z  x),  P2  =  (  j 2-?/ 2->-s 2)  an<^ 
P3  =  (£3,2/3,23)  such  that  Pj  ^  P2  #  P3,  these  points  determine  a  unique  plane 
in  space  through  the  determinant  equation: 


x  -  x1     y  -  yl     z  -  zj 

x2  ~   xl     Y2  ~   Yl     z2   ~    zl 

x2  -  xi  y~2  -  yi  z2  -  zi 


=  0 


To  simplify  the  equation,  we  replace  the  constant  differences  by  the  expressions: 


Clx  -  x2  — 
C2X  =  x3  - 
Cly  =  x2  - 
C2y  =  x3  - 
CI.  =  x9  - 


C2.  = 


Xo   - 
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Evaluating  the  determinant  using  the  diagonal  approach  yields: 

[(x  -  x1)ClyC2,  -  (x  -  x,)ClzC2y]  + 
[(y  -  yi)ClzC2x  -  (y  -  yJCl»C2j  + 
[(z  -  Zl)ClxC2y  -  (z  -  Bl)ClyC2j  =  0 

Solving  the  equations  for  x,y  and  z  in  terms  of  the  constant  expressions: 

A  =  ClyC2z  -  ClzC2y 
B  =  ClzC2x  -  C1XC2Z 
C  =    ClxC2y  -  ClyC2x 

D  =    -[AXl  +  Byi  +  CzJ 

The  Normal  to  the  Plane 

Once  the  parameters  A,B  and  C  have  been  determined,  the  solution  of  the 
linear  equation  for  any  normal  to  the  plane  is  straightforward.  Using  the  plane 
parameters  in  the  parametric  equation  for  the  normal  line  to  the  plane  and  using 
any  known  point  on  the  plane  (xkwn,ykwn.zkwn)  (the  midpoint  of  the  fractal 
triangle  for  example)  determines  a  normal  line  as: 

X  =  *kwn  +  c  A 
Y  =  ykwn  +  c  B 

Z   =    zkwn   +    c  C 

where  c  is  a  parameter  such  that  c  is  an  element  of  R.    By  varying  the  parameter 
c  we  can  solve  for  unique  points  on  the  normal  line  to  the  plane. 
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